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Quasi-Analytic  Model  of  OTHR  Clutter  from  Equatorial  Bubbles  in  the  Ionosphere 

Paul  A.  Bernhardt 
Plasma  Physics  Division 
Naval  Research  Laboratory 
Washington,  DC  20375 

Abstract 

Ground  clutter  for  the  over-the-horizon  radar  (OTHR)  sky  wave  becomes  Doppler  shifted 
because  the  ionosphere  through  which  the  radio  rays  are  propagating  changes.  One  source  of 
these  changes  is  ionospheric  bubbles  which  rise  vertically  through  a  horizontally  stratified 
plasma  near  the  equator.  The  rising  plume  structures  are  formed  when  gravity,  neutral  winds 
or  external  electric  fields  act  on  the  F-region  plasma.  The  effect  of  these  ionospheric 
disturbances  can  be  simulated  by  tracing  rays  through  physics  based  models  of  the  equatorial 
bubbles.  For  the  physics  based  models,  exact  solutions  for  internal  electric  potentials  are 
derived  assuming  linear  or  circular  symmetry  to  the  density  structures  imbedded  in  the 
background  plasma.  A  wide  variety  of  analytic  solutions  for  electric  potentials  are  found  for 
both  density  cavities  and  density  enhancements.  Quasi-analytic  solutions  for  the  transport  of 
the  bubbles  are  derived  using  the  continuity  equation  for  the  plasma  with  production  and  loss 
terms  neglected.  The  analytic  models  of  the  electric  fields  produce  incompressible  motion 
that  transports  the  locations  of  “plasma  cells”  but  do  not  change  the  density  of  the  plasma  in 
each  cell.  This  Lagrangean  approach  employs  a  time  dependent  coordinate  mapping  of  the 
undisturbed  layer  grid.  Using  internal  electric  potentials  of  the  bubbles  and  external 
polarizations  of  the  F-layer  as  a  whole,  a  transport  model  yields  tilted  plasma  plumes  that 
move  through  the  F-Region.  This  time-dependent  computer  model  provides  useful  plasma 
densities  in  a  fraction  of  the  time  for  fully  numerical  simulations.  The  electric  potential 
derived  in  the  models  can  be  directly  applied  to  the  ray  trace  computations  to  yield 
predictions  for  Doppler  shifts  in  the  unstable  ionosphere. 

I.  Introduction 

High  frequency  (HF)  over  the  horizon  radar  (OTHR)  systems  can  track  aircraft,  missiles  and 
ships  for  thousands  of  kilometers.  These  distances  are  achieved  with  sky-wave  propagation 
that  involves  reflection  and  refraction  in  the  earth’s  ionosphere.  Under  disturbed  ionospheric 
conditions,  the  radar  return  spectra  can  be  broadened  by  the  motion  of  structures  in  the 
ionosphere  where  the  radar  path  propagates.  Radar  signals  that  propagate  near  the  equator 
can  encounter  plasma  bubbles  that  are  formed  after  sunset  at  the  bottomside  of  the  F-layer. 
When  the  ionosphere  is  no  longer  illuminated  by  the  sun,  ion-electron  recombination  causes 
the  bottom  of  the  F-region  to  become  steeper.  Simultaneous  with  this  steepening,  the  layer  is 
lived  by  electric  fields  driven  by  the  dusk  terminator  neutral  winds.  This  process  causes  the 
F-layer  to  become  unstable  as  the  Rayleigh-Taylor  instability  forms  bubbles  that  form  on  the 
bottom  and  subsequently  rise  through  the  layer.  The  purpose  of  this  report  is  to  formulate  a 
model  of  these  bubbles  and  to  use  ray-trace  techniques  to  determine  the  effects  that 
equatorial  structures  with  have  on  OTHR  clutter. 
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OTHR  clutter  is  primarily  the  result  of  backscatter  from  the  ground  that  reflects  the  radar 
signal  back  to  the  source  with  a  range  delay.  Motion  of  the  ground  reflection  surface  such  as 
ocean  waves  can  induce  Doppler  spread  in  the  clutter  signal.  Even  if  the  ground  surface  is 
fixed,  spread  clutter  can  be  produced  if  the  ionosphere  changes  in  time  along  the  radar  path. 
The  sources  of  short  term  ionospheric  fluctuations  that  can  produce  Doppler  spread  in  the 
radar  clutter  are  (a)  decay  of  the  layer  after  sunset,  (b)  motion  of  the  layer  by  tidal  electric 
fields,  (c)  traveling  ionospheric  disturbances  caused  by  acoustic  gravity  waves  propagating  in 
the  neutral  atmosphere,  and  (d)  plasma  instabilities  that  cause  large  scale  irregularities  in  the 
plasma.  The  equatorial  bubble  is  the  most  import  irregularity  that  produces  clutter  at  low 
latitudes.  The  ionospheric  bubble  OTHR  clutter  study  proceeds  in  two  phases.  First  a  quasi- 
analytic  model  of  equatorial  bubbles  is  generated  that  (1)  replicates  all  the  features  of 
naturally  occurring  ionospheric  density  structures  and  (2)  can  be  used  to  calculate  radar  paths 
using  ray  tracing  through  the  electron  densities.  Rays  will  be  traced  from  a  ground 
transmitter  site  through  this  model  ionosphere  to  eventually  arrive  at  on  the  ground  for 
backscatter  along  the  same  path  to  the  transmission  point.  Each  ray  is  characterized  by  phase 
and  ground  path  lengths.  Temporal  changes  in  the  phase  path  produce  Doppler  shifts  in  the 
return  echo.  The  group  path  gives  the  round  trip  time  delay  that  will  be  assigned  to  a  range 
bin.  The  frequency  shift  in  each  range  bin  is  determined  from  the  time  rate  of  change  of  the 
phase  path  along  the  ray. 

II.  Overview  of  Ionospheric  Bubble  Models 

The  F-Region  ionosphere  can  become  unstable  if  a  density  perturbation  becomes  electrically 
polarized  by  external  forces  from  electric  fields,  neutral  winds,  and  gravitational  acceleration. 
Near  the  geomagnetic  equator,  gravity  can  act  on  the  plasma  attached  to  the  nearly  horizontal 
magnetic  field  lines  to  produce  unstable  conditions.  After  sunset  when  the  layer  is  lifted  by 
ambient  electric  fields,  the  bottom-side  steepens  and  plasma  bubbles  are  formed.  These 
bubbles  rise  through  the  layer  in  response  to  a  Rayleigh-Taylor  type  instability.  Also,  winds 
or  electric  fields  induce  electric  fields  in  both  density  cavities  and  enhancements  that  cause 
distortions  in  the  density  structures.  These  distorted  plasma  structures  are  responsible  for 
distortion  of  radio  propagation  which  lead  to  navigation  errors  and  outages,  communications 
systems  failures,  radar  clutter,  and  degradation  of  surveillance  systems.  The  modeling  of 
ionospheric  bubbles  or  density  enhancements  uses  computer  simulations  the  calculate  the 
effects  of  self  generated  electric  fields  (E)  that  are  driven  by  gravity,  neutral  winds  and 
external  electric  fields.  The  equations  for  these  simulations  can  be  solved  numerically  using 
(1)  flux  corrected  transport  algorithms  for  transport  of  plasma  and  (2)  direct  or  iterative 
solvers  of  the  non-separable  potential  equations  that  describe  the  self-generated  electric 
fields.  The  computational  time  for  solving  for  the  disturbed  ionosphere  is  often  prohibitive 
so  analytic  solutions  to  both  the  transport  and  potential  equations  are  useful.  Exact  analytic 
solutions  can  be  used  for  (1)  testing  of  the  numerical  algorithms  to  determine  errors  produced 
by  boundary  conditions  and  numerical  round-off  and  (2)  time-dependent  simulations  of 
typical  electron  densities  used  for  testing  tomographic  reconstructions  of  the  disturbed  F- 
layer  and  for  tracing  of  radio  ray  paths  through  the  region.  The  analytic  solutions  also  yield 
insight  into  the  conditions  for  production  of  ionospheric  bubbles. 
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The  computational  and  analytical  techniques  for  simulations  of  equatorial  bubbles  are 
compared  in  Figure  1.  Usually  numerical  models  proceed  as  illustrated  by  the  block  diagram 
given  in  Figure  la.  A  stratified  model  of  the  F-layer  is  perturbed  by  a  small  density 
disturbance.  Gravity  is  allowed  to  setup  an  electric  potential  in  this  plasma.  The  electric 
potential  is  obtained  with  a  numerical  solution  of  a  non-separable  elliptic  equation  using  a 
direct  solver  such  as  described  in  Appendix  A.  Once  the  electric  potential  is  obtained,  the 
plasma  is  transported  in  response  to  the  electric  fields  for  a  small  time  step.  A  non- 
dissipative  flux  corrected  transport  algorithm  (Zalesak,  Ossakow  and  Chaturvedi,  1982)  is 
then  used  to  incrementally  move  the  plasma  disturbance.  The  process  is  repeated  with  the 
generation  of  a  revised  electric  potential  followed  by  more  incremental  plasma  transport.  All 
of  these  processes  are  numerically  intensive  and  can  require  several  hours  of  computation. 
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Figure  1.  Block  diagram  of  (a)  numerical  and  (b)  quasi-analytic  algorithms  for  ionospheric 
bubble  modeling.  Both  approaches  use  an  equivalent  set  of  equations  but  apply  different 
solution  techniques  and  different  frames  of  reference. 
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The  quasi-analytic  model  for  the  equatorial  bubble  is  derived  using  the  three  steps  in  Figure 
lb.  This  model  is  computationally  faster  than  the  fully  numerical  model  and  is  more  suitable 
for  raytrace  analysis.  The  components  of  the  model  are  described  separately.  First  the 
electric  potential  is  defined  by  an  analytic  function  and  self-consistent  expressions  for 
electron  density  (or  Pedersen  conductivity)  structures  are  obtained  in  the  presence  of  electric 
fields,  neutral  winds  and  gravity.  Second,  the  plasma  transport  is  determined  using 
incompressible  motion  from  the  induced  electric  fields.  The  plasma  transport  is  derived  with 
the  analytic  electric  potential  distorting  the  coordinates  without  changing  the  density  in  each 
coordinate  cell.  To  realistically  model  ionospheric  bubbles,  the  effects  of  neutral  winds  as 
well  as  gravity  must  be  included.  The  polarization  of  the  layer  by  a  neutral  tilts  the 
ionospheric  plume  from  the  vertical.  The  third  step  is  to  adjust  the  parameters  in  the  analytic 
models  until  the  analytic  solution  given  in  step  1  matches  the  quasi-analytic  results  from  step 
2.  The  details  of  the  electric  potential  and  the  plasma  transport  processes  are  described 
Sections  HI  and  IV,  respectively.  Introducing  tilts  and  bifurcations  in  the  model  bubbles  is 
described  in  Section  V  and,  finally.  Section  VI  describes  normalization  of  the  model  to  yield 
realistic  time  and  spatial  scales.  Finally  in  Section  VII,  the  propagation  theory  for 
determining  Doppler  shifts  from  ionospheric  motion  is  described.  A  simple  example  for  the 
velocity  shifts  from  a  rising  ionosphere  is  used  to  test  the  model.  Future  work  will  use  the 
computer  models  derived  here  to  simulate  the  spread  of  the  OTH  radar  clutter  signals  by  both 
vertical  and  horizontal  motion  of  bubbles  in  the  equatorial  ionosphere. 

III.  Analytic  Models  for  the  Electric  Potential  in  a  Disturbed  Ionosphere 

The  equatorial  ionosphere  is  commonly  thought  of  as  a  uniform  layer  with  the  occasional 
imbedded  structure  or  bubble.  The  modeling  of  ionospheric  bubbles  uses  computer 
simulations  that  calculate  the  effects  of  self  generated  electric  fields  (E)  that  are  driven  by 
gravity,  neutral  winds  and  external  electric  fields.  The  equations  for  these  simulations  can  be 
found  in  a  number  of  papers  including  Bernhardt  [1988].  For  the  analytic  solutions 
considered  here,  the  background  ionosphere  will  be  uniform  in  the  horizontal,  x-  and  z- 
directions.  The  ambient  magnetic  field,  B,  is  aligned  with  the  z-axis.  The  altitude  variations 
of  the  undisturbed  ionosphere  will  be  represented  by  the  function  ne0(y)  where  y  is  the 
vertical  coordinate. 


The  layer  becomes  distorted  when  a  small  perturbation  grows  as  electric  fields  provide 
incompressible  perpendicular  motion  at  F-region  altitudes.  These  internal  electric  fields 
move  plasma  across  magnetic  field  lines  with  the  velocity 


v  = 


ExB 


VOxB 

Bo 


(1) 


where  v  is  the  velocity  perpendicular  to  the  ambient  magnetic  field  B  and  the  electric  field 
E  =  -VO  can  be  represented  as  the  gradient  of  a  scalar  electrostatic  potential  <b.  Other 
perpendicular  components  of  velocity  driven  by  pressure  gradients,  neutral  winds,  and 
gravity  can  be  neglected  because  the  plasma  in  the  F-region  ionosphere  is  magnetized.  This 
means  that  the  electron  and  ion  gyro  frequencies  are  much  larger  than  the  corresponding 
collision  frequencies  with  the  background  neutral  gas. 
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An  analytic  approach  is  derived  to  solve  for  the  electric  potential  for  localized  electron 
density  perturbations  driven  by  external  vector  fields.  Kelley  [1989]  gives  the  electric 
current  in  the  F-region  as 

J  =  [E  +  E0  +  (U+— )xB]ap  (2) 

Vi„ 

where  ap  =  6  n<;Vin  is  the  Pedersen  conductivity,  n^x,  y,  z)  is  the  electron  density,  Vjn  is  the 

B0 

ion  neutral  collisions  frequency,  Eo  is  the  external  electric  field  perpendicular  to  B, 
U  =  Uxx  +  Uyy  gives  the  normal  components  of  the  neutral  wind  vector,  g  =  -g0y  is  the 

gravitational  acceleration,  and  B  =  B0z  is  the  magnetic  field  vector.  Assume  that  the  normal 

electric  fields  are  constant  along  the  magnetic  field  in  the  z-direction  and  that  there  is  no  z- 
component  of  current  J.  For  specific  gravitational  accelerations,  neutral  winds  and  external 
electric  fields,  the  potential  equation  is  found  from  (2)  using 

|^+Vx.J  =  0  (3) 

dz 

Substitution  of  (2)  into  (3)  and  integration  along  the  z-direction  leaves  an  equation  for  the 
potential  in  the  perpendicular  (±)  x-  and  y-  directions. 

Vx  •  (ZpVx<D)  =  ZpVx2<l> + Vx0>  •  VxZp  =  Vx  •  j[E0  +  (U+— )xB]  apdz  (4) 

^in 

where  Ep  =  J*apdz  is  the  field-line-integrated  Pedersen  conductivity  and  E  =  — VX<I>  is  the 

o 

induced  electric  field.  For  simplicity,  the  driving  fields  E0,  U,  and—  are  assumed  constant 
in  space  and  time,  then  the  potential  equation  simplifies  to 


Vx2®  =  {[E0  +  (U  +  -^) xB] -  Vxd>}  •  Vx ln(E  )  =  {ET  -  Vx<t»}  •  Vx  ln(Ep)  (5) 

Vin 


g 

where  the  equivalent  electric  field  vector  defined  by  ET  =  E0  +  (U  +  —  )  x  B .  Given  a 

spatial  distribution  for  the  Pedersen  conductivity  (or  electron  density),  the  potential  is  usually 
obtained  numerically  from  the  non-separable  elliptic  equation  (5).  Often  iterative  solvers 
requiring  relatively  long  solution  times  or  direct  solvers  requiring  large  memories  are 
required  to  compute  this  solution. 

A  computational  alternate  approach  assumes  that  the  potential  is  given  and  (5)  is  used  to  find 
the  associated  electron  density.  For  this  solution,  only  Pedersen  currents  in  the  horizontal,  x- 
direction  will  be  considered  so 
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where  Etx  represents  the  equivalent  driving  fields  from  a  static  electric  field  in  the  positive  x- 
direction,  a  neutral  wind  in  the  positive  y-direction  and  gravitational  acceleration  along  the 
positive  y-axis.  The  sign  of  Etx  is  negative  for  the  downward  acceleration  of  gravity.  In  our 
notation,  the  growth  rate  for  the  Rayleigh-Taylor  instability  is  y  =  (-ETx/B0)/LN  where  Ln  is 

the  scale  length  of  the  gradient  on  the  bottom  side  of  the  ionosphere  [Zalesak,  Ossakow,  and 
Chaturvedi,  1982;  Sultan,  1996]. 

In  normalized  Cartesian  coordinates,  (4)  becomes 


a2<|>  3Log(zp)  a<!>  a2o  ,  9Log(zp)  a<|)  3Log(zp)  _ 


ax2 


ax 


ax  +  + 


dy  dy  dy 


ax 


=  0 


(7) 


where  often  gravitational  forcing  is  the  sole  contribution  to  Etx,  6  = 


<D 


^Tx  r0 


is  the 


dimensionless,  normalized  potential,  (X,  y)  = 


'x  y^ 
v.ro  ’  Toj 


are  the  normalized  coordinates,  and  ro 


is  a  constant  scale  factor  for  all  distances.  Note  that  (7)  has  many  self-similar  elements. 
Multiplying  the  X-  and  y-  coordinates  as  well  as  6  by  a  constant  scale  factor  (i.e.,  ro)  does 
not  change  the  equation.  Multiplying  Zp  by  a  constant  C0  also  yields  a  solution. 
Consequently,  if  <I>(x,  y)  and  Tp(X,y)  satisfies  (7)  then  so  do  the  pairs  of  functions 

^  X  V  XV 

r0  4>( — ,— )  and  C0  Tp( — ,— ).  Normalized  coordinates  (X,  y)  will  be  used  to  simplify  the 


ro  ro 


ro  ro 


notation  for  the  analytic  solutions. 


The  existence  of  analytic  solutions  for  (7)  was  discovered  by  examining  numerical  solutions. 
A  numerical  algorithm  for  non-separable  elliptic  equations  similar  to  (7)  was  written  using  a 
block  tri-diagonal  solver  of  the  algebraic  equation  derived  from  finite  difference 
approximations  to  the  partial  derivatives  (Appendix  A).  When  a  circularly  symmetric 

function,  Zp(r)  where  r  =  *Jx2  +  y2  ,  was  used  for  the  Pedersen  conductivity  it  was  found 

that  the  integral  of  the  resulting  electric  potential  along  X  was  also  circularly  symmetric.  In 
mathematical  terms, 


£P(r)  =>  J<I>(x,y)  dX  =  F,(r) 


(8) 


This  immediately  shows  that  the  form  of  the  potential  is  the  X  coordinate  multiplied  by  a 
circularly  symmetric  function  since 


(9) 


a>(x,y)  = 


F,(f)  F,(r)  dr  _  F,(r)  x  _ 

3x  dr  dx  dr  r  aU 


Functions  of  single  variables  such  as  r  can  be  solved  analytically. 

General  analytic  solutions  to  (7)  can  be  derived  by  making  simplifying  assumptions  about  the 
form  of  <f>  and  Log(Zp).  Numerical  simulations  for  symmetric  perturbations  in  x  with  x- 
directed  fields  given  by  (6)  yield  electric  potentials  that  are  odd  functions  of  x  with  the  form 

0(x,  y)  =  (a0  +  a,  x)  G[q(x,  y)]  (10) 


where  ao  and  a.\  are  constants,  q(x,y)  =  -y/x2+sy2  is  a  single  variable  representing  an 
elliptical  perturbation  for  the  potential  and  “s”  determines  the  polarization  of  the  coordinate 
ellipse. 

The  Pedersen  conductivity  (or  electron  density)  takes  the  form  of  an  elliptically  shaped 
perturbation  modulating  the  background  conductivity. 

Logics,  y)]  =  Lp[q(5E,y),s]+ loglS^Cy)]  (1 1) 


where  Lp(q,s)  is  the  natural  log  of  the  conductivity/density  perturbation 

is  the  integrated  conductivity  associated  with  the 


and  £ 


po(y)=  J- 


e  ne0(y,  z)vm 


r0dz 


'b 

horizontally  stratified  plasma  layer.  The  last  term  of  the  left  side  of  (7)  is  x-direction 
gradient  of  the  log-Pedersen  conductivity  which  drives  the  solution  for  the  potential. 
Consequently,  the  x  variation  through  the  function  q  is  required  to  obtain  useful  solutions  for 
the  potential. 


Substitution  of  (10)  and  (11)  into  (7)  and  solving  for  the  derivative  of  Lp(q)  yields  the 
equation 


q2G'(q)  +  s  +  sy8l°gl?po(y)]}  +  q3G'(q)  +  y 2  (s  -  l)s[qG'(q)  -  G'(q)] 


L'p(q,s)  =  - 


a0  +a,x 


3y 


xq2[a,G(q)  1]  +  q3Q,(q)  +  ?  2  (g  _  i)sqG'(q)] 


a0  +a,x 


(12) 

where  the  prime  (0  denotes  the  derivative  with  respect  to  q.  If  the  functions  of  y  vanish, 
(12)  may  be  integrated  directly.  The  y2  terms  vanish  only  if  s  =  0,  or  1. 
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General  solutions  for  an  extended  vertical  plume  imposed  on  a  horizontally  stratified 
ionosphere  are  considered  first.  For  s  =  0  the  potential  has  no  variations  in  the  altitude 
coordinate  (y )  and  the  solutions  for  (10)  and  (1 1)  are  given  by 


q  =  |x| 

4><0)(x,y)  =  (a0  +a,x)  G(q) 

L> ,  0x  _  2  a,  x  G'(q)  +  q  (a0  +  a,x)  G*(q) 

2[a,G(q)-l]+q(a0+aix)G#(«D 

2  a!  xG'(q)  +  q(a0+a,x)G'(q) 


E'())  (x,  y)  =  C0£  o(y)  exp 


-I 


x[a,  G(q)-l]  +  q(a0+a,x)G/(q) 


dq 


(13) 


where  Co  is  a  constant  chosen  to  give  the  background  conductivity  as  x  — » °° .  Exact 
solutions  for  pairs  of  Pedersen  conductivity  and  electric  potentials  from  (13)  are  easily  found. 
Table  I  gives  several  examples  these  pairs. 


Table  I.  Electron  Density  Disturbances  and  Compani 

ion  1-D  Potential  Functions 

Pedersen  Conductivity  Function,  E(p0)  (x,  y) 

Electric  Potential,  <P(x,  y) 

£po(y)exp(bx2) 

(a0  +a,  x)exp(-b  x2) 

exp(bx2)  +  2b(a0  +a,  x)x-a, 

2po(y)(i+|x|b)2 

aj  x 

(l+|x|b)2  +  a,(b-l)|x|b  -a, 

i+ 1  x  r 

Sp0(y)2cosh[bx]2 

(a0  +a,x)  sech[b  x] 

1  -  2a,cosh[b  x]  +  cosh[2  b  x]  +  2  b(a0  +  a,x)  sinh[b  x] 

Figure  2.  One-Dimensional  Pedersen  conductivity  associated  with  the  1-D  analytic  electric 
potential  using  the  parameters  ai  =  -0.5  and  b  =  3.  The  topology  of  the  solution  remains 
unaffected  by  the  choice  of  model  parameters. 
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An  example  of  the  second  function  in  Table  I  is  illustrated  in  Figure  2  with  the  parameters  ai 
=  -0.5  and  b  =  3.  The  conductivity  trough  (Figure  2a)  centered  along  the  y-axis  has  a  ridges 
that  increase  in  amplitude  as  the  parameter  “ai”  is  increased.  If  ai  >  0,  the  trough  is  replaced 
by  a  Pedersen  conductivity  enhancement  as  the  sense  of  the  potential  (Figure  2b)  is  reversed. 
The  parameter  “b”  simultaneously  controls  the  steepness  of  the  walls  on  the  conductivity 
irregularity  and  the  spatial  decay  of  the  potential  function.  With  separation  of  variables,  even 
more  general  solutions  to  (7)  can  be  found  in  Cartesian  coordinates  (Appendix  B). 

General  Solutions  for  Circular  Holes  in  an  Ionosphere  with  simple  vertical  structure  is 
considered  next.  For  s  =1  the  density  disturbance  is  circularly  symmetric  with  radius  r 
around  the  ( x ,  y )  origin.  The  solutions  from  (12)  require  ao  =  0  and  ai  =  1  with  the  result 


r  =7x2  +y2 
0(1)(x,y)  =  xG[r] 


(3  +  y 


di°gPpo(y)] 


Lp  (?,!)  =  ■ 


dy 


-}G'(r)+rG'(r) 


G(r)  +  rG'(f)-l 


(14) 


To  eliminate  any  dependence  on  of  y  in  (14),  the  background  Pederson  conductivity  takes 
the  functional  form 


Ep0(y)  =  C0ym  where  m  and  Co  are  a  constants.  With  this  substitution,  (14)  becomes 


Lp(f,l)  =  - 


(3  +  m)G'(f)  +  r  G*(r) 
G(r)  +  ?  G'(r)  - 1 


(15) 


which  is  identical  to  (13)  with  ao  =  0,  ai  =  1  and  m  =  -1.  The  corresponding  formula  for  the 
spatial  variation  of  Pedersen  conductivity  is 


Kl)(z>y)=c 


0  ym  exp| 


-I- 


(3  +  m)G/(f)  +  fG<'(r) 
G(r)  +  rG'(r)-l 


(16) 


In  Appendix  C  using  separation  of  variables,  all  of  the  solutions  for  (7)  in  cylindrical 
geometry  are  derived. 


The  one-dimensional  (s=0)  and  two-dimensional,  circularly-symmetric  (s=l)  expressions  are 
similar.  With  ao  -  0  and  ai  =  1,  the  rational  polynomial  function  of  the  form 


G(q)  = 


l  +  qb 


used  in  (10),  and  the  corresponding  potentials  from  (13)  and  (14)  are 
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(17a) 


O(0)(x)  = 


x  a 

i+itf 


and  <i>(1)(x,y)  = 


x  a 


l  +  (x2  +  y2)b/2 


(17b) 


where  a  and  b  are  constants.  Analytic  solutions  can  be  obtained  from  (13)  through  (16). 


Derivatives  of  the  density  perturbations  become 

Tv  abqb-'[l  +  b  +  (l-b)qb] 

L (q,0)  = - r3— 1 - c -  .  ,  where  q=  x 

(1+  qb)[a-a(b-l)qb-(l  +  qb)2] 


t  i\  &  b  r  [2  +  b  +  m  +  (2-b  +  m)r  ]  ,  ^  [Z 5  ~ 

L(r,l)  = - h- - £ - ,  where  r  =  Jx 2+y2 

v  '  (1+  rb)[a-a(b-l)rb-(l  +  rb)2]  v 

The  corresponding  Pederson  conductivity  expressions  are 

(l+|x|b)2 


i?(M)-i5(S)-i_a  ^ 


^1)(x,y)=  exp 


a(l  +  m) 
i - : - 


l  A, 


2tan' 


B,  -2rt 

A, 


-  sign(a)  n 


c0  ym  (i  +  fb)2 

1  -  a  -  B,r b  +  r2b 


(18a) 

(18b) 


(19a) 

(19b) 


where  A,  =  yj- a[4b  +  a(b  - 1)2]  and  B,=-a(b-l)-2  and  r  =  -y/x2  +  y2  .  The  constants  of 

integration  are  chosen  to  yield  the  background  density  at  large  distances  where  x  — »  «> .  The 
physically  acceptable  solutions  have  b  >  0. 

Two  types  of  conductivity  structures,  cavities  and  enhancements,  are  described  by  (19a  and 


-4b 


=  aMin  >  a  >  0 ,  the  plasma  structure  is  a  cavity  centered 


19b).  In  the  parameter  range  2 

at  x  =  0.  These  limits  are  found  by  solving  for  A\  =  0.  As  parameter  “a”  approaches  the 
value  of  “amin”  the  sides  of  the  density  cavity  becomes  steeper.  With  a  =  ami,,,  the  wall  of  the 


cavity  is  located  at  radius  is  r  = 


^  fb  +  m  +  2 
r  = 


b  +  1 
b-1 


For  b  >  m  +  2,  the  cavity  has  a  ridge  located  at 


b-m-2 


With  0  >  a  >  1,  a  conductivity  enhancement  is  found  at  the  origin.  This  enhancement  can 
represent  the  increased  Pedersen  conductivity  produced  by  an  artificial  ion  cloud  from 
Barium  or  similar  material  released  in  the  sunlit  ionosphere.  As  a  — »  1 ,  the  sides  of  the 
density  enhancement  become  steeper.  Solutions  exist  for  all  values  of  b  >  0  but  if  b  >  1,  the 
potential  vanishes  at  large  distances.  If  a  >  1,  then  the  solutions  in  (19a  and  19b)  become 
complex  and  are  not  physically  possible.  The  maximum  upward  velocity  for  the  potential 
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<f>(1)(x,y)  in  (17b)  is  V'  =  -  — ^  ^  =  -  ^Tx  .  If  a  >  1,  then  the  conductivity  would 

B0  3x  B0 

move  with  a  velocity  larger  than  the  Etx/Bo  velocity  of  the  driving  force,  which  is  not 
possible.  For  instance  if  only  a  vertical  wind  Uy  is  considered  in  the  equation  (6)  for  Etx, 
then  the  upward  velocity  is  VyO  =  a  Uy.  An  unreasonable  value  of  a  >  1  would  permit  the 
conductivity  enhancement  to  rise  faster  than  the  neutral  wind  driver. 


The  restrictions  on  the  ranges  for  the  potential  amplitude  “a”  indicate  that  not  all  electric 
potentials  correspond  to  a  physical  density  or  Pedersen  conductivity  structure.  For  a  given 
force  on  the  plasma  from  external  electric  fields,  neutral  winds  or  gravity,  the  induced 
potential  is  determined  by  the  gradients  on  the  wall  of  density  cavity  or  enhancement.  These 
gradients  are  physically  limited  by  infinite  steepness  and  the  amplitude  of  the  potential  is  a 
maximum  at  this  limit.  Thus,  for  the  solution  to  (4),  a  given  physical  density  structure  will 
always  correspond  to  a  potential  function.  The  magnitude  of  a  potential  function  can  be 
increased  to  the  point  that  there  is  no  corresponding  plasma  density  function. 

The  one  dimensional  expression  (19a)  can  represent  a  horizontal  density  modulation  that 
uniformly  changes  the  background  density  of  a  stratified  ionosphere.  These  may  be 
produced  by  horizontally  traveling  acoustic-gravity  waves  can  act  as  seeds  for  equatorial 
bubbles.  The  elongated  shapes  of  these  modulations  are  illustrated  by  the  example  in  Figure 
2.  The  elongations  can  be  found  in  nature  as  the  extensions  of  an  ionospheric  plume  below 
its  top.  The  horizontal  electric  field  vectors  calculated  as  gradients  of  the  potential  yield 
vertical  plasma  transport.  This  transport  is  normal  to  the  density  gradients  and,  consequently, 
no  net  change  in  the  densities  is  produced.  The  electric  fields  near  the  top  of  the  bubble  are 
the  primary  drivers  for  plasma  transport. 

The  two-dimensional  solutions  to  the  potential  equation  are  more  useful  than  the  one¬ 
dimensional  solutions.  The  expression  for  ^(x,  y)  in  (19b)  describes  a  plasma  disturbance 
with  two-dimensional  structure.  The  conductivity  (or  electron  density)  vanishes  in 
E(p°(x,y)  at  y=0  unless  m=0.  Figure  3  illustrates  three  examples  of  the  analytic  density 

cavities  and  the  associated  electric  potential  for  a  uniform  background  using  m  =  0.  By 
changing  the  parameters  in  the  analytic  model,  a  wide  variety  of  density  structures  is 
obtained. 


The  background  plasma  variation  can  be  approximated  using  nonzero  values  of  m.  The 
Pedersen  conductivity  from  (19b)  vanishes  at  the  y  =  0  boundary  if  m  >  0.  With  Ep  =  0  at 


the  lower  y=0  boundary,  the  potential  equation  (7)  reduces  to 


3£p  30 

3y  3y 


=  0 . 


To  satisfy  this 


condition,  dO/dy  =  0  because  vertical  gradient  3£p/3y  is  nonzero-positive  at  the  lower 
boundary.  This  condition  is  automatically  built  into  the  analytic  expression  for  the  electric 
potential  because  of  the  y 2  symmetry  of  d>(x,y)  =  x  G(^Jx2  +  y2) . 
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Pedersen  Conductivity 


Normalized  Potential 


b 


a  =  0.5 


Figure  3.  Analytic  results  for  density  cavities  and  enhancements  in  a  uniform  background 
(i.e.,  m=0).  The  densities  and  potentials  are  computed  using  (a)  b  =  2 ,  a  =  0.99  a  Min  =  - 
5.92,  and  (b)  b  =89,  a  =  0.5  a  Min  -  -0.48,  (c)  b  =  4,  a  =  0.5.  The  changes  in  the  parameters 
yield  either  (a)  a  cavity  with  steep  sides,  (b)  a  ridge  around  the  cavity  or  (c)  a  peaked 
enhancement. 
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Two  quantitatively  different  ionospheres  can  become  polarized  with  the  same  potential 
variation.  Figure  4  illustrates  two  solutions  of  (19b)  with  identical  parameters  for  a  and  b  but 
with  (a)  m  =  1  for  a  linear  background  profile,  and  (b)  m  =  0.25  for  a  forth-root  of  y  profile. 
As  seen  by  the  solutions  in  (19a  and  19b),  a  family  of  Pedersen  conductivity  structures  can 
be  associated  with  a  single  electric  potential  (or  field)  distribution.  This  non-uniqueness 
property  can  be  exploited  for  modeling  the  evolution  of  the  density  structures 


Pedersen  Conductivity  Normalized  Potential 


a  =  -1.78 
b  =  4 


Figure  4.  Density  cavities  imbedded  in  a  linear  and  inverse-quartic  variation  for  the  vertical 
profile  of  the  plasma  conductivity.  Both  structures  yield  the  same  electric  potential.  The 
conductivity  and  the  y-directed  electric  field  goes  to  zero  a  the  lower  boundary. 
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One  use  of  the  analytic  solutions  to  the  potential  equation  (7)  is  to  provide  test  cases  for  the 
numerical  solutions  of  the  same  equation.  Numerical  solutions  are  influenced  by  (1)  the 
finite  difference  approximations  to  the  derivatives,  (2)  the  boundary  conditions,  and  (3)  the 
convergence  error  for  iterative  solvers.  A  direct  numerical  solver  was  used  to  solve  for  the 
potential  (Appendix  A).  The  Pedersen  conductivity  was  specified  using  the  second,  two- 
dimensional  analytic  model  in  (19)  with  the  parameters  a  =  -0.53,  b  =  4,  and  m  =  1.  The 

^  x  a, 

analytic  potential  0(1,(x,y)  = - - - r^from  (17b)  is  compared  in  Figure  5  with  the 

1  +  (x  +  y  ) 

numerical  solutions  for  the  potential  with  a  variety  of  boundary  conditions. 

For  uniform  boundary  densities  as  illustrated  in  Figure  3,  doubly  periodic,  doubly  fixed  (or 
Dirichlet)  [with  the  boundary  potential  set  to  0  =  0],  or  doubly  derivative  (or  Neumann) 
[with  90/dx  =  0  and  dO/dy  =  0  normal  to  each  boundary]  specifications  work  well.  The 
non-uniformities  for  the  Pedersen  conductivities  at  the  boundaries  in  Figures  4  and  5  require 
care  in  the  selection  of  the  boundary  conditions  for  numerical  solutions.  The  Pedersen 
conductivity  and  corresponding  analytic  potential  are  shown  in  Figures  5a  and  5b, 
respectively.  As  illustrated  by  Figures  5c  and  5d,  inaccurate  solutions  are  obtained  for  the 
potential  with  boundaries  specified  as  doubly  periodic  and  doubly  fixed/Dirichlet  (0  =  0), 
respectively.  These  two  boundary  conditions  force  equal  potentials  at  the  top  and  bottom 
where  the  real  densities  are  different.  Neumann  (i.e.,  zero  derivative)  boundary  conditions 
yield  a  useful  solution  (Figure  5e).  Another  accurate  solution  is  obtained  by  using  a 
Neumann  (i.e.,  zero  derivative)  boundary  at  the  bottom,  a  Dirichlet  (i.e.,  zero  potential) 
boundary  at  the  top  and  periodic  boundaries  at  the  sides  of  the  solution  space  (Figure  5f). 
Table  II  lists  the  maximum  potential  for  each  solution  using  a  64  x  32  grid.  All  of  the 
numerical  solutions  with  the  correct  shape  (Figures  4e  and  4f)  yield  a  computed  potential  that 
is  about  1 1  %  less  than  the  actual  values.  This  reduction  is  the  result  of  the  finite  difference 
approximations  for  the  derivatives  in  (7). 

As  the  number  of  mesh  points  is  increased,  the  numerical  solution  becomes  more  accurate. 
Table  III  shows  the  effect  of  the  grid  size  on  the  maxima  of  the  computed  potentials.  The 
numerical  solutions  use  an  nx  by  ny  grid  in  the  x-  and  y -directions,  respectively.  In  all 
cases,  the  doubly  derivative  or  Neumann  boundaries  yield  slightly  better  solutions  than  the 
mixed  boundary  solution  with  doubly  periodic  in  the  x -direction  and  Neumann/Dirichlet 
(derivative/fixed)  boundary  values  in  the  y -direction.  These  examples  illustrate  that  the 
analytic  solution  pair  (17b)  and  (19b)  provides  an  easy  way  for  testing  (1)  the  utility  of  the 
numerical  solutions  with  various  boundary  conditions  and  mesh  sizes  and  (2)  that  the  error  in 
the  numerical  solution  vanishes  as  the  grid  becomes  denser. 
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Pedersen  Conductivity  Analytic  Potential  Solution 


Figure  5.  Comparison  of  analytic  and  numerical  solutions  for  the  electric  potential  with  a 
density  cavity  at  the  bottom  side  of  the  stratified  plasma.  The  potential  is  set  up  with  either 
downward  gravity  and  winds  or  horizontal  background  electric  fields.  The  magnetic  field  is 
in  the  z-direction  perpendicular  to  the  horizontal  x-axis  and  vertical  y-axis.  Analytic 
solutions  for  (a)  the  normalized  Pedersen  conductivity  and  (b)  the  normalized  analytic 
potential  are  compared  to  the  numerical  solutions  for  the  potential  with  (c)  x-periodic  and  y- 
per iodic  boundaries ,  (d)  x-fixed/Dirichlet  and  y-fixed/Dirichlet  boundaries,  (3)  x- 

derivative/Neumann  and  y-derivative/Neumann  boundaries,  (f)  x-  derivative/Neumann  at 
bottom,  x-fixed/Dirichlet  at  top  and  y-periodic  boundaries.  The  latter  two  boundary 
conditions  yield  numerical  solutions  that  approximate  the  analytic  model. 
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The  analytic  solutions  provide  a  quantitative  test  of  numerical  algorithms  for  the  non- 
separable  PDE’s  that  describe  electric  potentials  in  the  ionosphere.  All  numerical  solvers  use 
the  approximations  of  (1)  finitely  spaced  solution  samples  (2)  boundaries  at  finite  distances. 
The  analytic  solution  provides  infinite  spatial  resolution  with  no  boundary  limits.  Test  cases 
can  be  set  up  to  show  the  effects  on  the  numerical  solutions  of  boundary  placement  and 
boundary  types  as  well  as  discrete  sampling  of  the  coordinate  system. 

Table  II.  Error  Analysis  for  Numerically  Derived  Potentials 


Solution 

x-Boundary 

Condition 

y-Boundary 

Condition 

Maximum 

Potential 

Maximum 

Error 

Shape 

Fig.  3 
Panel 

Analytic 

— 

— 

0.304 

0 

Correct 

b 

Numerical 

Periodic 

Periodic 

1.167 

0.863 

c 

Numerical 

Fixed 

Fixed 

0.154 

0.150 

Incorrect 

d 

Numerical 

3<I>/3x  =  0 

3O/3y  =  0 

0.273 

0.031 

Correct 

e 

Numerical 

30/3x  =  0 

3<b/9y  =  0 /Fixed 

0.271 

0.033 

Correct 

— 

Numerical 

Fixed 

9<I>/9y  =  0 /Fixed 

0.269 

0.035 

Correct 

— 

Numerical 

Periodic 

3<X>/9y  =  0 /Fixed 

0.269 

0.035 

Correct 

f 

Table  III.  Mesh  Size  Affects  on  Numerically  Derived  Potentials 


Number  of  Cells 
in  x -Direction: 
nx 

Number  of  Cells 
in  y -Direction: 

% 

Maximum  Potential  (%  Error) 

Derivative 

Boundaries 

Mixed 

Boundaries 

32 

15 

JEBSHBSEM 

64 

32 

EmnEESH 

128 

64 

EEEERim 

ISEU30SM 1 

256 

128 

I'KMfclMI 

In  summary,  exact  analytic  solutions  have  been  found  for  the  nonlinear  potential  equation 
commonly  used  for  determination  of  electric  fields  for  ionospheric  plasma  irregularities. 
These  solutions  provide  easy  means  to  calculate  the  distributions  of  plasma  conductivity 
associated  with  analytic  models  for  the  electric  potential.  These  solutions  can  represent 
plasma  depletions  or  enhancements  depending  on  the  model  parameters.  The  analytic 
examples  demonstrate  that  a  large  family  of  density  structures  can  correspond  to  identical 
electric  potentials.  Also,  if  the  amplitude  of  a  potential  is  too  large,  there  will  not  be  any 
corresponding  electron  density  structure.  One  use  of  the  analytic  solutions  is  to  test 
numerical  techniques  to  solve  for  the  electric  potential  associated  with  arbitrary  distributions 
of  electron  density  under  the  influence  of  gravity,  winds,  and  ambient  electric  fields.  The 
next  step  to  complete  the  quasi-analytic  bubble  model  is  to  use  the  exact  forms  of  the  electric 
potential  derived  in  this  section  to  provide  transport  of  the  plasma  structures.  This  is 
described  in  the  next  section. 
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IV.  Quasi-Analytic  Solutions  for  Plasma  Transport  in  a  Disturbed  Ionosphere 

Plasma  densities  evolve  by  transport,  compression,  production  and  loss.  These  processes  are 
contained  in  the  continuity  equation 

^  +  V(nv)  =  P-L  (20) 

where  n  is  density,  P  is  production  and  L  is  loss.  Numerical  or  analytic  solutions  to  (20)  can 
be  based  on  either  the  Eulerian  or  the  Lagrangean  form  of  the  equations  [Richtmyer  and 
Morton,  1967].  These  forms  are  equivalent  except  the  Lagrangean  form  describes  where 
each  bit  of  fluid  cam  from  originally.  Equation  (20)  is  the  Eulerian  from  of  the  equation  of 
continuity.  The  Lagrangean  from  is  derived  for  the  special  case  of  incompressible  flow. 

The  compressibility  of  a  plasma  is  given  by  the  term  V  •  v .  When  this  term  is  zero,  the 
plasma  is  termed  incompressible.  Using  (1)  the  compressibility  of  the  F-region  plasma  is 
found  to  vanish  because 


V  •  v  =  -V  •  (V<D  x  B)/B02  =  [-B  •  (V  x  Vd>)  +  VO  •  (V  x  B)]/Bj  =0  (21) 

where  a  constant  B  is  assumed. 


Expanding  (20)  with  (21)  yields  the  compressionless  form  of  the  continuity  equation 


an  _  fa 

—  +  v-Vn  = 

at 


\ 


— +  v-V 

vdt  j 


Dn  T 
ins —  =  P-L 
Dt 


(22) 


where  Dn/Dt  is  called  the  total  derivative.  The  total  derivative  of  the  density  moves  with  a 
small  volume  element  (Ax,  Ay,  Az)  in  the  velocity  field.  During  this  process  a  plasma 

element  at  (x,  y,  z)  is  mapped  to  another  location  (x’,  y’,  z’)  in  an  increment  of  time  At .  The 
incremental  mapping  equation  is  given  by 

x'  =  x  +  v(x)At  (23) 


where  the  vectors  have  components  x  =  (x,  y,  z)  and  v  =  (vx,  vy,  vz).  The  derivative  form  of 
(23)  simplifies  (22)  so  that 


dn(x,  t) 
dt 
dx(t) 


=  P(x,t)-L(x,t) 


(24) 


dt 


=  v(x,t) 
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If  production  and  loss  are  neglected,  then  the  electric  fields  simply  move  volume  elements  of 
plasma  in  space  but  the  density  in  each  element  remains  unchanged.  Equations  (24)  are  the 
Lagrangean  form  of  the  continuity  equation. 

An  analytic  formulation  is  developed  to  describe  the  equatorial  bubbles  in  terms  of  a 
mapping  function  that  distorts  the  ionospheric  layer  according  to  the  second  equation  in  (24). 
The  mapping  function  is  usually  determined  with  a  numerical  simulation  that  calculates  the 
electrostatic  E  fields  as  a  function  of  time  and  space  as  an  ionospheric  bubble  or  irregularity 
is  formed  in  the  F-layer.  Substitution  of  these  fields  into  (1)  yields  the  transport  velocities 
and  the  second  equation  in  (24)  can  be  solved  to  provide  the  motion  of  the  density 
coordinates.  The  transformation  of  coordinates  by  this  process  is  given  by 

x(t)  =  M[x(t0),t-t0]  and  x(t0)  =  M_'[x(t),t-t0]  (25) 

and,  neglecting  production  and  loss,  the  electron  densities  are  given  by 

n(x,t)  =  n(M_1[x,t-t0],t0)  (26) 


where  the  map  M"1  transforms  the  distorted  coordinates  back  to  the  initial  coordinate 
locations. 

Taking  the  magnetic  field  to  be  aligned  with  the  z-direction,  the  electron  density  fluctuations 
are  assumed  to  vary  in  the  2-dimensional  coordinate  system  (x,  y).  The  coordinate  transform 
map  is  given  as 


x(t)' 

/— \ 
o 

1 

+-» 

© 

o 

& 

X 

S 

l _ 

_y(t)_ 

1 

© 

VS 

o 

r+ 

1 

f-t* 

O 

where  y  is  altitude,  x  is  zonal  distance  in  a  flat  earth  system  and  (xo,  yo)  are  the  initial  values 
for  these  coordinates  at  time  to. 

The  mapping  function  M(x,  y,  t)  must  be  one-to-one  and  invertible  and  be  the  identity  map 
where  the  induced  electric  potential  is  zero  and  the  plasma  densities  are  unchanged. 
Substitution  of  (1)  into  (24)  yields 

dl(0=_V®xB 

dt  Bq 

Assuming  uniformity  along  B  in  the  z-direction,  the  differential  equations  governing  the 
coordinated  transformations  are 


ax(x0,y0,t)  _ _ 1  9Q(x,y,t)  ^  ay(x0,y0,t)_  1  3d>(x,y,t) 

3t  B0  dy  9t  B0  dx 


(29) 


18 


Differentiating  (29)  by  xO  and  yO  and  using  the  Poisson’s  equation 


VE  =  V(-VO)  =  -(^+^) 
ox0  oy0 


yields  the  equations  for  the  velocities 

d(dx/9t)  _  8^/9t)  and  3(3x/dt)  8(9y/9t)  _ 

dx0  dy0  dy0  dx0 


(30) 


(31) 


The  areas  between  curves  of  constant  xO  and  yO  are  preserved  in  the  transformation  and  the 
trajectories  of  the  coordinate  transformation  follow  contours  of  constant  <3>(x,y).  The  electric 
potential  is  setup  by  gravity,  neutral  winds  and  external  electric  fields. 


Consider  the  coordinate  transformation  provided  by  the  analytic  potential  from  (17b)  in  the 
previous  section 


a>(x,y) 


(x/ro)  ax  ib  ETx 
l  +  [(x/r0)2  +  (y/r0)2]b/2 


*  ax  r0  Et* 

l  +  fb 


(32) 


where  coordinates  x  and  y  are  the  Cartesian  coordinates  normalized  by  a  constant  scale 
factor  ro  and  constants  ax  and  b  define  the  shape  of  the  potential. 

The  plasma  transport  velocities  are 


9x  _  b  x  y  r  h™2 
9t  "  (l  +  fb)2 
ay  1-bx2  fb-2  +  fb 

at-  (l  +  rb)2 


(33) 


where  t  =  &x  Etx  t  and  r  =  Jx2+y2  . 
roBo 

The  maximum  upward  velocity  at  the  (x  =0,  y=  0)  origin  is  given  by  =  —  ^Tx , 

B0 

independent  of  the  potential  shape  parameter  “b”.  Note  that  in  our  example  of  the  downward 
gravity  vector  driving  the  transport,  the  parameter  Etx  from  (6)  is  less  than  zero  and  a  value 
of  ax  <  0  is  required  to  yield  an  upward  velocity.  The  parameter  ax<0  denotes  a  density 
cavity  and  consequently  the  center  of  the  cavity  is  expected  to  rise  against  gravity.  If  the 
parameter  “ax”  were  greater  than  zero,  the  center  of  the  density  enhancement  would  fall  as 
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expected  under  the  influence  of  gravity.  For  the  rest  of  the  discussion,  only  density  cavities 
will  be  considered. 


The  sides  of  a  cavity  fall  under  the  influence  of  gravity.  The  minimum  downward  velocity  of 


v  -  axEtx 

vy>  - 


(b-1)2 


=  -V. 


(b-1)2 


B  4b  *  4b 

coordinate  system  is  distorted  by  the  potential  inside  a  conductivity  cavity  so  that  the  center 
cells  move  upward  and  the  side  cells  move  downward. 


is  found  at  x  =  ± 


'b+lV* 


b-1 


and  y=0.  The  Cartesian 


The  analytic  model  potential  with  b  =  4  was  previously  illustrated  in  Figure  3b.  The 
corresponding  vector  field  for  the  plasma  velocities  (Figure  6a)  shows  the  central  uplift  of  the 
plasma.  As  a  result  of  this  flow,  after  normalized  time  t  =  1  the  initially  square  cells  become 
mapped  according  to  the  results  shown  in  Figure  6b.  The  horizontal  (red)  and  vertical  (blue) 
grid  lines  become  distorted  by  the  vortex  flow  from  the  potential.  Note  that  the  area  in  each 
plasma  cell  remains  constant  during  this  process. 


a> 

13 

c 


T3 

L. 

8 

o 


Normalized  Velocity  Field  Distortion  of  Square  Grid  Cells 
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Normalized  x-Coordinate 
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Figure  6.  Computed  (a)  plasma  velocities  and  (b)  mapping  of  Cartesian  coordinates  by  the 
analytic  model  for  the  electric  potential  shown  in  Figure  3b.  The  longest  velocity  vector  has 
a  E 

a  magnitude  - —and  the  spatial  coordinates  are  normalized  by  the  scale  length  ro .  The 

B 

coordinate  map  is  illustrated  after  the  flow  velocities  flow  have  operated  on  the  plasma  for 

i 

or  when  normalized  time  7  =  1 . 


time  t  = 


f  r  V1 
aE* 

\Bro 


The  primary  feature  of  ionospheric  bubbles  is  that  they  rise  leaving  an  elongated  cavity  of 
reduced  plasma  density  often  referred  to  as  a  “plume”.  As  the  plume  rises,  it  carries  the 
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along  electric  potential.  Using  the  upward  velocity  for  the  potential  in  (32)  of  =  — — — 

B0 

the  analytic  model  for  the  potential  becomes 
0(x,  y,  t)  =  0>(x,  y  -  V^O)  =  - 


xax  ETx 


1  + 


'x*  ( 


vroy 


y  V 


V 


b/2 


\vr0  r0  J 


(34) 


With  this  potential  formula,  the  plasma  drift  velocities  are  found  from  (29)  to  be 


b-2 


dx  ^  b  [x2  +(y- 1)2]  2 

=x(y-t)-^  - - ^-r 


l  +  [x2+(y-t)2]z 


b-2 


9y  l  +  [(l-b)  xa+(y- t)2][x2+(y-t)a1  ■ 

3T  j,+[^+(y-V]nJ 


(35) 


where  the  time  has  been  normalized  by  t  =  t  V^/rg  and  [x(0),  y(0)]  =  [x0,y0]  are  the  initial 

conditions  at  time  t  =  t0  =  0 .  This  equation  is  the  base  for  the  quasi-analytic  transport  for 

production  of  ionospheric  bubbles.  By  using  this  equation,  it  is  assumed  that  the  electric 
potential  shape  is  totally  specified  with  fixed  parameters  Etx,  a,  b,  and  ro.  The  only  temporal 
variation  of  die  potential  is  that  it  rises  with  a  constant  speed  given  by 


The  potential  in  the  trailing  portion  of  the  plume  is  neglected  because  any  vertical  flow 
below  the  bubble  only  operates  on  horizontal  gradients  and,  consequently,  there  is  no  net 
plasma  transport  within  this  region. 

In  the  reference  frame  of  the  rising  potential,  yp  =  y  - 1  and  the  differential  equations 
become 


dx 

8T 


=  xyp 


b-2 


b  [x2  +y,l  2 


and 


<8 

dt 


p  _ 


=  -l  + 


b-2 


1+[(1- 

b)  x2  +y,][ 

x2  +Yp* 

2 

1 

|l  +  [x2+y2 

f 

(36) 
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The  coordinate  transformation  if  found  by  integrating  (36)  for  (x,yp)  with  an  initial  value  of 
(x0,y0)  and  then  finding  y  =  yp  - 1  .  Further  simplification  is  obtained  by  transforming  to 
spherical  coordinates  where  x  =  rp  Cos0p  and  yp  =  rp  Sin  0P  .  With  this  substitution,  (36) 
becomes 


5?i 


p  _ 


at 

30P 

at 


rp  sin  0P 
l  +  fp 

Cb-l 


rp'  cos0p(l+b  +  rp ) 


(37) 


(l  +  rp  )  where  the  initial  conditions  at  t  =  0  are  rp(0)  =  ->jxl  +  yl  and 


0P(O)  =  tan-1(y0/x0)  . 


The  coordinate  transform  starts  with  an  electric  potential  at  altitude  y  =  t0  =  0  and  lets  this 

potential  distort  the  coordinates  until  time  t,  when  the  potential  is  at  altitude  y  =  t, .  The 
starting  and  stopping  times  and  altitudes  are  critical  in  defining  the  coordinate  distortions. 
For  describing  equatorial  bubbles,  these  starting  altitude  must  be  transformed  to  the  location 
where  the  bubbles  starts  to  form  on  the  bottom  side  of  the  ionosphere  and  the  stopping 
altitude  yields  with  the  location  of  the  bubble  at  time  t,  =  t,  .  This  renormalization  is 
described  later. 

The  properties  of  this  coordinate  transform  can  be  examined  at  the  center  where  x  =  0  and 
(35)  simplifies  to 


=  0  with  x(0)  =  x(t)  =  0 
8t 


8y 

at 


— - - with  y(0)  =  y0 

i+[(y-t)2F 


(38) 


The  solutions  of  (38)  are  found  from  the  nonlinear  equation 

ySign®o)+^lfl=|y0|+lMl 

1-b  1-b 


(39) 


The  limiting  forms  for  the  asymptotic  solutions  to  (39)  are 
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y  =  y0for  y0<-l 
y  =  y0  +  tif  t  <  tc  1 

y  =  y0  +  tcif  t  >  tcJ 

y  =  y0  for  0  <  y0  and 


for  - 1  <  y0  <  0  and  tc 

t<y0 


y  =  t  for  y0  =  0 


(-y0)‘~b 

b-1 


y=t  + 


(b-i) 


l/b 


for  1  <  y0  and  t  =  y0 


y  =  t  + 


(40) 
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Figure  7.  Mapping  of  the  (a)  y-axis  and  (b)  a  square  coordinate  grid  with  a  rising  potential 
that  starts  at  the  origin.  Using  equation  (34)  with  b  —  4,  the  coordinates  above  the  initial 
center  of  the  potential  are  swept  into  a  narrow  band  that  rises  with  the  upward  bubble 
velocity.  The  snapshot  at  7  =  6  shows  transport  of  the  ionosphere  upward  at  the  center  and 
downward  at  the  sides. 
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This  maps  of  the  coordinate  distortions  are  illustrated  in  Figure  7  using  the  potential  with 
index  b  =  4.  The  temporal  variations  in  the  normalized  altitude  are  given  in  Figure  7a.  The 
spatial  mapping  of  the  grid  at  time  t  =  6  is  shown  in  Figure  7b.  The  ionosphere  is  nearly 
undisturbed  below  the  starting  altitude  (y  =  0)  of  the  rising  electric  potential.  The 

ionosphere  is  also  undisturbed  for  times  ( t  <  y0 )  when  the  potential  is  below  the  undisturbed 
ionosphere.  As  the  potential  rises,  it  sweeps  up  all  the  coordinates  and  carries  them  in  a 
narrow  layer  at  just  above  the  center  of  the  potential  function  where  y  =  t .  This  type  of 
coordinate  transformation  can  be  used  to  simulate  the  rising  equatorial  bubble.  At  time  t,  the 
potential  rises  to  be  centered  at  an  altitude  y,  =  t, . 


Figure  8.  Coordinate  compression  at  the  top  of  the  plasma  bubble.  For  the  potential 
parameters  b  =  4,  the  normalized  y  -coordinates  are  mapped  to  a  thin  shell  about  0.4  to  0.45 

above  the  center  of  the  rising  potential  function  located  at  y  =  t. 


The  horizontal  coordinates  near  the  center  of  the  rising  potential  map  to  a  thin  shell  centered 
at  the  altitude  y  =  t .  An  approximation  to  this  map  is  given  by  the  last  equation  in  (22).  The 
evolution  of  this  coordinate  “compression”  is  illustrated  in  Figure  8  for  normalized  times 
between  t  =  1  and  20  using  the  parameter  b=4  for  the  electric  potential.  The  offset  (  y  - 1 ) 
from  the  center  of  the  rising  bubble  increases  monotonically  but  slowly  as  the  initial 
coordinate  covers  a  much  larger  range. 
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The  normalized  density  gradient  at  the  edge  of  the  bubble  is  ^  so  is  the 

dy  9y0  dy  dy 

coordinate  compression  factor.  Analytically,  this  factor  is  given  by 


dyp  _  (  y0Yi+yb 

dy  L  y  J  1+yo 


(41) 


This  factor  is  plotted  in  Figure  9  with  b  =  4  for  a  wide  range  of  times  and  initial  coordinate 
altitudes.  The  compression  factor  increase  with  time  easily  attaining  values  greater  than  10 
or  100.  With  this  compression,  a  bottomside  ionospheric  gradient  becomes  greatly  amplified 
as  the  bubble  rises. 


Figure  9.  Density  gradient  enhancement  factor  for  the  coordinate  compression  at  the  top  of 
the  rising  bubble  with  index  b  -  4.  The  compression  becomes  enhanced  with  increases  in  the 
normalized  time  ? .  Outside  the  effect  of  the  potential  this  factor  is  unity.  This  processes 
yields  the  steepened  gradients  at  the  sides  of  the  ionospheric  bubble. 
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To  obtain  the  electron  density  or  Pedersen  conductivity  at  time  ( t ),  cell  coordinates  must  be 
obtained  for  the  original  cell  that  was  transported  to  the  current  position.  The  model 
ionospheric  bubble  is  therefore  formed  by  operating  on  a  horizontally  stratified  layer  with  the 
inverse  of  the  coordinate  transformation  defined  by  (35).  This  inverse  may  be  obtained  with 
by  (1)  interpolation  of  the  solution  to  (35)  or  (2)  by  reversing  time  for  the  solution  of  (35). 
Figure  7  shows  the  location  of  the  coordinates  at  the  current  time  for  each  normalized 


starting  position  (x0,  y0)  = 


y  o 


The  interpolation  process  is  hampered  by  the 


v  ro  r«; 

distortion  of  the  coordinate  density  cells.  In  the  mapping  shown  by  Figure  7b,  the  square  cell 
with  comers  at  (x0,y0)  =  (0,  2),  (0.1,2),  (0.1,  2.1),  and  (0,  2.1)  are  mapped  to  the  elongated 


region  with  comers  (0,  6.421),  (0.1075,4.374),  (0.1112,  4.443),  (0,  6.424)  at  t  =6.  With 
linear  interpolation,  the  point  (x,  y)  =  (0.055,  5.4)  inside  the  elongated  region  is  determined 
to  map  to  the  initial  position  (x0,y0)  =  {0.0505,  2.0364).  This  is  in  error;  the  correct 
location  for  this  position  is  (x0,  y0)  =  (0.3329,  -0.1784).  Beside  being  prone  to  error,  this 
technique  requires  interpolation  on  a  non-uniform  mesh  or  numerical  solution  for  the  inverse 
of  interpolation  on  an  a  uniformly  spaced  mesh  of  the  coordinate  locations  at  time  t  =  0 .  For 
good  accuracy  and  ease  of  computation,  interpolation  or  numerical  inverse  solutions  should 
be  avoided. 


Since  the  electric  potential  does  not  evolve  with  time,  the  inverse  coordinate  mapping  is 
easily  achieved  with  time  reversal.  Consider  a  cell  with  location  (xj,  yi)  at  time  t  =  fi. 
Running  time  backwards  to  time  t  =  0  yields  its  starting  point  (xo,  yo).  If  the  parameters  Vyo 
and  b  are  constant,  the  time  reversal  solution  is  most  easily  obtained  by  replacing  t  with  —  t 
in  (35)  through  (37).  This  process  yields  the  map  represented  by  (26). 

The  equations  for  the  inverse  coordinate  map  are 


dgo 

3t 


=  -x0(y0  +  t) 


b 

Xo2 +(yo +  02] 

b-2 

2 

I'* 

[x02+(y0  +  t)2 

H-[(l-b)x;+(y,  +  t);fe+(y,  +  y1  ■ 


(42) 


The  center  of  the  potential  function  starts  at  altitude  y0  =  t,  and  then  (42)  solved  as  the  center 
of  the  potential  falls  to  an  altitude  y0  =  0 .  For  this  reason,  the  inverse  coordinate  map 
equations  are  integrated  from  t  =  - 1,  to  t  =  0 .  The  initial  conditions  for  (42)  are 

x0(-t,)  =  xlandy0(-t,)  =  y1  (43) 
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Normalized  Distance,  ? 


The  inverse  coordinate  map  is  illustrated  in  Figure  10  for  the  initial  grid  and  the  distorted 
inverse  grid  at  several  times.  This  map  is  used  to  determine  the  origin  of  a  coordinate  cell 
and  the  initial  electron  density  or  Pedersen  conductivity  in  that  cell.  As  an  example  of  using 
this  inverse  map,  the  mapping  of  specific  point  (0.055,  5.4)  is  directly  obtained  with  (42)  to 
yield  the  correct  initial  location  (0.3329,  -0.1784). 


Figure  10.  Inverse  coordinate  map  showing  source  regions  for  density  gradients. 

No  analytic  solution  to  (42)  could  be  found  but  the  numerical  solution  to  the  coupled 
ordinary  differential  equations  is  fast  and  efficient  using  standard  explicit  methods.  The 
Lagrangean  inverse-map  is  consequently  called  quasi-analytic  because  no  finite-difference 
approximations  are  applied  to  the  spatial  (x,  y)  coordinates  but  discrete  steps  in  the  time 
dependent  solutions  of  (42)  is  required.  Applying  the  inverse  map  to  any  analytical  model  of 
a  uniform  ionosphere  can  efficiently  provide  the  density  at  any  point  in  space  and  time  that 
can  be  obtained  without  the  use  of  interpolation  or  the  numerical  solution  of  two- 
dimensional  partial  differential  equations. 
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As  an  example,  the  inverse  coordinate  map  is  applied  to  an  F-layer  described  by  the 
following  formula  for  a  modified  Chapman  layer 


Ne(y)  =  Ne0  Exp[  1  -  z  -  Exp(z)] 


z  = 


y-HP 

H0 


(44) 


H0  -  H0I  +  [0.5  +  Tan-‘(^-^)  /  n]U02 

“i 

The  parameters  Neo,  Elp,  Ho,  Hoi,  H02  and  Hi  control  the  shape  of  the  layer.  The  analytic 
simulation  uses  peak  density  Neo  =  106  cm"3,  peak  altitude  Hp  =  400  km,  bottom-side  scale 
height  Hoi  =  20  km,  top-side  scale  height  H02  =  50  km,  and  transition  scale  Hi  =  10  km.  This 
simple  layer  model  has  a  steep  bottomside  representative  of  the  equatorial  ionosphere.  The 
conversion  from  normalized  coordinates  is  y  =  y  r0  +  yc0  where  yco  is  the  starting  altitude  of 
the  electric  potential  at  time  t  =  0. 
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Figure  11.  Simulation  of  a  bubble  formed  in  the  ionosphere  using  a  simple  formulation  for 
the  electric  potential.  The  parameters  for  the  model  are  b=4,  and  ro  =  30  km. 
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Applying  the  inverse  coordinate  transformation  (42)  to  the  ionospheric  profile  yields  the 
bubble  evolution  illustrated  in  Figure  1 1 .  The  coordinate  distances  are  determined  using  a 
scaling  factor  ro  =  30  km  and  the  potential  function  rises  through  the  layer  starting  at  yc0  = 
370  km  altitude.  The  potential  function  index  is  arbitrarily  set  to  b  =  4  for  this  example.  The 
normalized  time  coordinate  t  =  t  V^/r0  =  t  ax  ETx/(Br0)  is  used  because  the  parameter  “ax 

has  yet  to  be  specified.  The  allowable  values  for  a  were  previously  given  after  (19)  as 
—  4b 

- -  <  a  <  0  .  With  a  larger  value  of  parameter  “ax  ETx”,  the  vertical  velocity  of  the 

(b  1) 

bubble  increases  and  the  absolute  time  (t)  in  the  simulation  is  reduced  for  a  fixed  normalized 
time.  Figure  1 1  illustrates  that  the  analytic  model  using  the  rising  potential  yields  a  quasi- 
analytic  solution  that  resembles  numerical  solutions  requiring  much  more  computation  time. 
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Figure  12.  Effect  of  electric  potential  index  b  on  the  model  ionospheric  bubble.  The 
normalized  time  for  each  solution  is  t  =  3.75 . 

Whereas  the  parameter  “ax”  controls  the  bubble  rise  rate,  the  “b”  determines  the  size  of  the 
bubble.  The  parameter  “b”  may  have  any  value  greater  than  unity.  Larger  values  of  b  yield 
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electric  potentials  with  larger  gradients  at  the  edges.  Increasing  b  increases  the  region 
affected  by  the  electric  potential.  Figure  12  illustrates  the  ionospheric  bubble  for  several 
values  of  this  index.  Values  of  the  potential  index  in  the  range  2  <  b  <  4  seem  to  yield 
reasonable  descriptions  of  the  rising  ionospheric  bubble. 

V.  Bubble  Tilts  and  Ambient  Shear  Flow 


Both  observations  and  numerical  computations  [Zalesak,  Ossakow  and  Chaturvedi,  1982] 
have  demonstrated  that  neutral  winds  can  tilt  the  ionospheric  plumes  off  vertical.  When 
acted  on  by  a  zonal  neutral  wind  Ux,  a  vertical  electric  field  Ey  is  produce  by  interactions 
with  the  background  electron  density  (or  Pedersen  conductivity).  Also,  a  perturbation 
electric  field  Ejy  is  produced  by  polarization  of  the  plasma  density  bubble  by  the  neutral 
wind.  These  two  processes  work  simultaneously  to  affect  the  tilt  the  bubbles.  First,  the 
background  neutral  wind  induces  large  scale  plasma  motions  in  the  zonal,  eastward  direction. 
The  shear  in  this  large-scale  plasma  drift  will  cause  the  regions  of  lower  background 
Pedersen  conductivity  to  lag  behind  the  regions  of  maximum  Pedersen  conductivity.  Second, 
polarization  of  the  ionospheric  bubbles  by  the  neutral  wind,  which  moves  faster  than  the  bulk 
motion  of  the  plasma,  will  produce  horizontal  motion  that  will  cause  the  density  depletions  to 
move  in  the  opposite  direction  of  the  neutral  wind  relative  to  the  bulk  plasma  drift.  This  is  a 
well  known  property  of  plasma  holes  as  they  respond  to  neutral  winds  [Bernhardt,  1988], 
Both  of  these  processes  are  easily  incorporated  in  the  quasi-analytic  bubble  model. 

Vertical  gradients  in  the  background  density  yield  vertical  gradients  in  the  induced  electric 
field  and  vertical  shears  in  the  horizontal  plasma  drifts  produced  by  these  fields.  This 
process  is  captured  in  (2)  assuming  that  the  vertical  currents  vanish  with  the  result 

Jy  =  [Ey  -  UXB0]  crp  =  0  (45) 


The  field  line  is  divided  into  the  F-region  where  the  wind  is  uniform  and  the  E-region  where 
the  neutral  wind  will  assume  to  vanish  [Zalesak,  Ossakow  and  Chaturvedi,  1982].  Calling 
the  integrated  Pedersen  conductivity  in  these  two  regions  EF  and  ZE  respectively,  the  vertical 
electric  field  profile  is  given  by 


E,(y) 


U,B„S„(y)  = 

EE+EF(y)  dy 


(46) 


where  Oo  is  the  polarization  potential  of  the  background  plasma. 
The  resulting  plasma  velocity  is  given  by  (24)  with  the  result 

—  =  — =  ux  sp(y)  =U  - — - =  V s 

3t  B0  xEE+ZF(y)  xf(y/rO)  +  l 


(47) 
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where  f(y/r0)  =  f(y)  = 

2F(y/r0) 


(48) 


This  is  the  horizontal  velocity  of  the  plasma  in  which  the  bubble  is  imbedded.  Usually  the 
wind  shear  variation  is  small  compared  to  the  average  bulk  motion. 

In  normalized  coordinates,  this  wind  shear  equation  for  the  background  horizontal  motion 
becomes 

dS-V.s_U,  -  u.  1  -y  m  (49. 

9t  V*  V„SE  +  LF(y/r0)  V„  1+Sy)  ,sW 

Integration  of  (49)  yields  the  simple  coordinate  map  from  this  large  scale  plasma  motion 


x(t)  =  jVxS(y)  dt  =x0  +  t  VxS(y) 


(50) 


where  xq  is  the  initial  coordinate  at  time  t  =  0 . 


Neutral  Density  (crrr3)  Electron  Density  (1 06  cnrr3) 


Figure  13.  Profiles  from  neutral  atmosphere  and  electron  density  models  used  for  the 
sample  computations  of  electric  polarization  that  tilts  the  equatorial  bubbles. 
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Conductivity  (m2/Ohm)  (m/s) 

Figure  14.  Model  profdes  of  field-line  integrated  Pedersen  conductivity,  uniform  neutral 
wind,  large-scale  horizontal  plasma  drifts,  and  horizontal  wind  relative  to  the  drifting 
plasma. 

If  the  E-region  Pedersen  conductivity  is  zero,  then  the  plasma  will  move  horizontally  with 
the  neutral  wind  speed  Ux.  A  finite  coupled  with  vertical  shears  in  the  F-region  Pedersen 
conductivity  gives  a  shear  structure  to  the  horizontal  plasma  motion.  The  tilts  of  the 
equatorial  plume  structures  will  be  computed  using  the  Lagrangean  approach  to  the  plasma 
transport  with  both  imbedded  potential  given  by  (34)  and  the  large  scale  distortions  described 
by  (50). 

As  an  example  of  this  wind  induced  shear  in  the  horizontal  plasma  drift,  a  uniform  neutral 
with  Ux  =  100  m/s  was  used  to  polarized  the  plasma  layer  given  by  (44)  with  a  peak  density 
of  106  cm'3.  The  model  neutral  atmosphere  and  electron  density  profiles  illustrated  in  Figure 
1 3  are  used  to  derive  the  equatorial  profile  of  the  field-line  integrated  Pedersen  conductivity 
shown  in  Figure  14.  A  dipole  magnetic  field  model  of  the  form 
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and  Be  = 


(51) 


_  2  H0  Cos9 

Dr  = - - - 


H0  Sin0 


is  used  for  the  magnetic  field  lines  where  H0  —0.311  10"5  Tesla.  Figure  15  illustrates  the 
distortion  of  a  square  mesh  using  the  coordinate  transform  defined  by  (47).  Application  of 
this  transform  after  the  transform  illustrated  in  Figure  7  will  tilt  the  bubble  to  the  left  (west) 
side  of  the  simulation.  The  simulation  for  Figure  15  used  a  fixed  E-region  conductivity  ZE 
that  was  one-tenth  the  maximum  value  of  the  F-layer  Pedersen  conductivity  Ep  . 


Figure  15.  Horizontal  coordinate  distortion  by  wind  induced  plasma  shear  flow  after  500 
seconds  of  plasma  motion. 


The  rising  bubble  will  be  caught  in  the  sheared  plasma  flow  to  provide  a  tilt  to  the  bubble. 
To  model  this  tilt,  the  quasi-analytic  transport  model  will  be  modified  assuming  that  the 
bubble  follows  a  trajectory  that  is  a  combination  of  the  vertical  rise  velocity  VyO  and  the 
horizontal  shear  velocity  given  by  (47).  The  equations  for  the  trajectory  of  the  center  (xc,  yc) 
of  the  bubble  is 

and  ^  =  V,s(yc)  (52) 

Assuming  that  the  bubble  starts  to  form  at  the  initial  location  (x,  y)  =  (0,  yco),  then  the 
solution  of  (52)  becomes 


yc  =  yco  +  VyO1  and  xc  =  Jvsx(yc0  +  VyoO  dfi  (53) 

0 
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In  normalized  coordinates  this  solution  for  the  center  trajectory  of  the  bubble  becomes 


yc  =  t  and  xc  =  jVxS(x)  df  (54) 

o 

where  distance  is  normalized  by  ro,  velocities  are  normalized  by  Vyo  and  time  is  normalized 
by  ro/Vyo  as  before. 

The  simple  coordinate  map  (47)  neglects  horizontal  flow  from  internal  polarization  of  the 
bubble.  This  flow  is  the  result  of  vertical  polarization  fields  generated  by  polarization  of  the 
bubble  structure  by  the  neutral  wind.  This  process  has  been  described  in  Section  II  except 
that,  along  with  vertical  gravitational  acceleration,  the  horizontal  neutral  wind  induces  an 
electric  potential  that  causes  horizontal  motion  of  the  bubble  center.  From  (6),  the  vector 
electric  field  from  these  forces  is 

E0  +  (U  +  — )xB  =  -Et  (56) 

Vin 

The  horizontal  neutral  wind  in  the  rest  frame  of  the  plasma  bubble  is  Ux  -  Vxo  where  Vxo  is 
the  horizontal  velocity  of  the  center  of  the  bubble.  The  electric  field  associated  with  this 
neutral  wind  in  this  frame  is 


ETyy  =  (Ux-Vx0)B0y 


(57) 


by  the  definition  in  (56).  In  the  rest  frame  of  the  plasma,  the  relative  velocity  of  the  electric 
potential  at  altitude  y  =  yc , 


vxR(yc) 


ETy  =  i  ao, 

b0  dy 


(58) 


where  <E>]  is  the  potential  from  internal  polarization. 

The  response  of  the  plasma  to  this  drive  field  is  dependent  on  the  physical  structure  of  the 
bubble.  For  the  tilted  bubble  model  with  internal  polarization,  a  single  parameter  scaling  of 
the  ambient  drifts  is  used.  The  horizontal  drift  velocity  of  the  bubble  center  is  assumed  to 
have  the  form 


Vx0(y)  =  (1  -  ay)  VxS  =  VxS  +  V*  with  V*  =  -  ay  VxS  (59) 

where  ay  is  a  constant  analogous  to  ax  used  in  the  previous  section  to  determine  the  bubble 

rise  rate  from  gravity.  The  parameter  ay  determines  the  bubble  velocity  in  a  rest  frame  of  the 
background  plasma.  For  density  depressions,  ay  is  less  than  zero. 

The  wind  and  plasma  drift  profiles  in  Figure  14  illustrate  this  polarization  effect.  Figure  14b 
shows  that,  near  the  peak,  the  relative  wind  in  the  plasma  rest  frame  is  about  10%  of  the  total 
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plasma  drift.  The  parameter  ay  controls  the  relative  velocity  of  the  bubble  in  the  background 
plasma.  If  ay  =  0,  the  bubble  drifts  with  the  background  plasma  as  if  there  were  on  internal 
polarization  of  the  bubble.  If  ay  =  -1,  the  background  horizontal  plasma  drift  adds  to  the 
wind-induced  relative  motion  produced  by  the  polarization  fields  inside  the  bubble  giving 
Vxo  =  2  Vxs  and  Vxr  =  -Vxs-  In  normalized  coordinates  (59)  becomes 


V  = 

vx0 


'  xO 


V 


=  (l-a,)VxS=VlS+V, 


1  xR 


yO 


(60) 


Lin  et  al.,  [2005]  measure  relative  zonal  propagation  of  equatorial  plasma  bubbles  that  is 
consistent  with  this  internal  polarization  effect.  Using  the  definition  of 


a  =  -V  N 

ay  vxR'  vxS 


(61) 


where  Vxr  is  the  measured  bubble  velocity  relative  to  the  background  plasma  drive  Vxs  ,  the 
observations  of  Lin  et  al.  [2005]  yield  time  dependent  values  of  ay  from  1.7  at  20  magnetic 
local  time  (MLT)  to  0.3  at  22  MLT.  Lin  et  al  [2005],  however,  attribute  this  effect  to  strong 
coupling  between  atmospheric  gravity  waves  and  the  Rayleigh-Taylor  instability.  This  paper 
asserts  that  polarization  potentials  of  the  density  depletions  set  up  by  the  relative  neutral  wind 
can  explain  the  enhanced  motion  of  the  bubbles. 

The  technique  for  bubble  modeling  is  based  on  the  motion  of  the  analytic  electric  potential 
along  a  trajectory.  The  tilted  bubble  rises  along  the  path  defined  by  the  velocity 
%  =  Vx0x  +  V^y .  With  both  background  drift  and  internal  plasma  motion,  the  dynamics  of 
the  center  of  the  plasma  bubble  potential  are  given  by 

^  =  V,,  and  ^  =  Vl0  =  (1  -  ay)V^  (62) 


with  the  initial  conditions  [xc(0),yc(0)]  =  [0,  yc0  ]  . 

For  the  tilted  bubble  model,  the  internal  horizontal  and  vertical  electric  fields  are  considered 
along  with  of  the  overall  motion  of  the  background  plasma.  With  the  internal  field 
assumption,  the  analytical  potential  function  becomes 


0(x,  y,  t)  =  O0  (y)  +  O,  (x  -  xc ,  y  -  yc0  -  V^t,  0) 

[x  -  xc(t)]  Vyp  -  [y  -  yc0  -  Vypt]  (yc0  +  V^t) 


=  ®0(y)+B0- 


1  + 


xc(t) 


\2 


y-YcO-Vyot 


ib/2 


(63) 


which  is  identical  to  (34)  with  <I>o,  Vxr  and  Xc  equal  to  zero.  The  induced  plasma  bubble 
velocity  in  the  horizontal  direction  is 
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(64) 


vx0  = 


J_3<D 
B0  3y 


V  +  V 

vxS  T  vxR 


at  the  position  x  =  xc(t)  and  y  =  yc0  +  V^t .  The  added  variables  tilt  the  electric  potential  off 

vertical  so  that  the  head  of  the  bubble  can  flow  against  the  ambient  drift  of  the  background 
plasma. 

The  coordinate  shift  xp  =  x-xc(t)  and  yp  =  y - yc0 - V^t  translates  the  potential  into  the 
reference  frame  of  the  bubble  center  with  the  result 


<D(x,y,t)  =  0)()(y)  +  B0 


XpVyO  -  yp  YrtCycO  +  Vyot) 

i+MJ’ 


(65) 


where  rp2=xj+yj  . 


The  potential  function  given  by  (65)  is  substituted  into  the  Lagrangean  transport  equations 

dx(x0,y0,t)  =  1  3<D(x,y,t)  ay(x0,y0>t)  =  1  9<B(x,y,t) 

3t  B0  dy  3t  B0  3x 

The  resulting  coordinate  mapping  equations  for  the  rising  potential  in  sheared  plasma  flow  is 
^  =  V.o(y.)  =  V<0(ydl+Vj0t) 

dx  _b  r^VyptXpVyo-yp  VxR(yc0  +  Vy0t)]  +  [l  +  rpV]VxR(yc0  +  Vy0t)  w 

U  b  -bV  +  vxsW  k00) 

*  (l  +  rpbr«7 

ay  V<,tl  +  r,V]-l>r„t-Vx1,[xcV'<,-y.Va,(ya,+V<,t)] 

*  (lnV? 

incorporating  both  internal  and  external  forces  on  the  bubble. 

In  normalized  coordinates,  the  transport  equations  become 
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(67) 


^  =  vx0(yc),  yc  =  t, 

ax  b  ^b-2y [x  -  VsR  (yc )  y]  +  [  1  +  ]  VsR  (y c )  ~ 

at  [i +?]2  18  w 

ay  l  +  fb-b^2x[x-VsR(yc)y] 

at  [i+^b]2 

where  once  again  distance  is  normalized  by  ro,  velocity  is  normalized  by  Vyo,  and  time  is 
normalized  by  ro/Vyo .  The  velocity  functions  are  all  related  to  the  normalized  plasma  shear 
by theparameter ay withVx0(y)/(l-ay)  =  -VxR(y)/ay  =  VxS(y)  . 

Before  computing  the  coordinate  mapping,  the  trajectory  of  the  center  of  the  potential 
function  (xc,  yc)  is  found  by  solving  the  first  ordinary  differential  equation  in  (66)  or  (67). 
For  this  calculation  and  all  the  follow,  the  shear  velocity  shown  in  Figure  15  was  used  in  (60) 
with  a  several  values  of  ay.  Figure  16  illustrates  the  model  results  for  the  motion  of  bubble 
center  as  a  function  of  the  parameter  ay.  With  ay  =  0  the  potential  will  rise  and  drift  east 
reaching  a  maximum  distance  on  the  topside  where  the  wind  induced  drift  vanishes.  As  the 
parameter  ay  is  increased  toward  unity,  the  internal  polarization  of  the  bubble  inhibits  its 
horizontal  motion  in  the  background  plasma  flow. 

To  illustrate  the  results  from  these  coordinate-mapping  equations,  the  parameter  ay  is  set  to 
one-tenth  so  the  potential  function  moves  horizontally  with  the  ambient  plasma  flow  at  each 
altitude  and  is  slightly  retarded  by  internal  polarization.  When  ay  =  0.1,  then  VxR(y)  = 
Vxs/10,  Vxo  -  9  Vxs/10.  The  initial  conditions  for  the  coordinate  map  are 
[xc(0), x(0), y(0)]  =  [0, x0,y0]  at  t=0.  All  of  the  results  are  displayed  in  normalized 

coordinates.  The  normalization  parameters  ro  =  30  km,  and  v^  =  100  m/s  can  be  used  to 
translate  the  solutions  back  to  physical  space.  The  bubble  moves  with  the  vertical  rise 
velocity  from  internal  electric  fields  set  up  by  gravity.  By  setting  the  internal  polarization 
parameter  ay  nearly  to  zero,  the  horizontal  motion  is  primarily  with  the  background  plasma 
but  there  is  a  small  retardation  from  the  internal  fields.  The  shear  function  was  put  into 
normalized  form  from  the  physics  coordinates  using  the  definition  y  =  (y  -  y^)/^  where  yco  - 
370  km.  With  the  parameter  b  =  4,  the  ordinary  differential  equations  given  in  (67)  are 
integrated  in  time  to  yield  the  quasi-analytic  solutions  for  the  Lagrangean  coordinate 
distortions  shown  in  Figure  17.  The  center  of  the  potential  function  as  derived  from  the  first 
equation  in  (67)  is  shown  by  the  green  curve  in  each  figure.  The  plume  structure  below  the 
top  of  the  bubble  becomes  caught  up  in  the  ambient  flow  to  from  a  backwards  “C”.  The 
successive  images  in  Figure  17  are  normalized  shown  for  normalized  times  of  2,  4,  and  6. 
Using  the  normalization  factor  r0/Vyo  =  300  seconds,  the  absolute  times  for  the  tilted  bubble 
coordinate  maps  are  (a)  600,  (b)  1200,  and  (c)  1800  seconds. 
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Figure  16.  Curves  of  the  trajectory  for  the  analytic  potential  function  used  to  create  the 
equatorial  bubble.  With  parameter  ay  =  0,  the  internal  polarization  of  the  bubble  is 
neglected  and  the  potential  function  is  transported  by  the  full  action  of  the  plasma  shear.  As 
ay  is  increased,  the  zonal  motion  of  the  center  of  the  bubble  potential  is  reduced.  With  ay  = 
1,  the  internal  polarization  at  the  head  of  the  bubble  completely  cancels  the  plasma  drift  of 
the  ambient  layer. 


38 


Normalized  Distance, 


o* 


1  2  3  4  5  1 

Normalized  Distance,  x 


Figure  17  Coordinate  map  for  distortions  from  a  rising  (b  =  4)  bubble  in  a  sheared 
background  plasma  flow.  The  green  line  shows  the  trajectory  of  the  center  of  the  bubble  with 
ay  -  0.1.  The  plume  becomes  curved  as  it  is  caught  in  the  ambient  plasma  flow. 


The  coordinate  transform  equations  in  (67)  yield  the  normalized,  destination  coordinates 
[x(t),y(t)]  from  the  given  initial  coordinates  [x0,y0].  To  determine  a  mapped  density  at  a 

given  location,  the  time-inverse  transformation  to  determine  the  initial  coordinates  for  a 
coordinate  cell  that  is  transported  to  a  given  location.  This  inverse  transform  has  already 
been  discussed  in  the  previous  section  for  bubbles  with  out  tilts.  This  inverse  map  are 
obtained  by  replacing  t  with  -t  in  (67).  The  center  of  the  potential  function  starts  at  location 
Xo  =  Xc(t,)  and  y0  =  t,. 


The  solution  for  the  inverse  coordinate  transform  proceeds  in  two  steps.  First,  the  initial 
equation  in  (67)  is  solved  to  yield  the  function  Xc(t)  for  the  horizontal  displacement  of  the 
bubble  potential  function.  Next,  the  system  given  by  (68)  is  integrated  as  the  potential 
follows  a  trajectory  [xc(t),  yc(t)]  to  end  at  the  origin  of  the  normalized  coordinate  system 

where  x0  =  0  and  y0  =  0 .  The  curves  in  Figure  16  show  the  [xc(t),  yc(t)]  trajectories  as  a 
function  of  the  internal  polarization  factor  ay.  As  previously  discussed  with  (42),  the  inverse 
coordinate  map  equations  are  integrated  from  t  =  —  t,  to  t  =  0 .  The  initial  conditions  for  (63) 
are 


x0(-t,)  =  x,  and  y0(-t1)  =  y, 


(69) 
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Figure  18  Computed  examples  of  (a)  inverse  coordinate  map  at  t  =  4  and  (b)  corresponding 
ionospheric  bubble  densities  at  t  -  1200  seconds  for  the  electric  potential  moving  through  a 
sheared  plasma  with  small  internal  polarization  (ay  —  0.1).  The  parameters  for  the 
simulation  are  identical  to  those  used  to  generate  Figure  11  which  can  be  used  for 
comparison  to  illustrate  the  effects  of  the  zonal  wind  on  tilting  the  bubble. 

The  inverse  coordinate  map  for  the  rising  bubble  the  sheared  plasma  flow  is  illustrated  in 
Figure  18a  in  normalized  coordinates  at  the  normalized  time  t=4.  Using  (68),  this 
coordinate  map  presents  the  origin  of  the  plasma  cells  that  have  been  transported  to  the  west 
by  the  wind  induced  horizontal  transport  and  the  gravity  induced  vertical  transport.  The 
parameter  a y  is  set  to  zero  so  the  horizontal  transport  by  internal  polarization  of  the  equatorial 
bubble  is  neglected. 

When  the  inverse  coordinate  map  is  applied  to  the  background  plasma  layer,  the  electron 
densities  are  to  form  the  tilted  plume  structure  in  Figure  18b.  The  depletion  at  the  center  of 
the  bubble  is  the  result  of  incompressible  convection  from  the  bottom  to  the  topside  of  the 
layer.  In  this  model,  the  reduced  density  channel  extends  over  100  km  down  though  the  layer 
to  the  bottom  of  the  F-region.  The  absolute  spatial  dimensions  are  derived  using  a  scale 
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length  ro  =  30  km  and  a  base  height  yco  =  370  km.  As  with  Figure  16,  the  time  normalization 
factor  is  300  seconds. 

VL  Parameter  Normalization  in  the  Analytic  Models  of  the  Disturbed  Ionosphere 

In  the  simulations  of  the  previous  three  sections,  the  electric  potential  was  fixed  as  it  rose 
through  the  plasma  layer  forming  bubble  structures.  The  next  step  in  the  quasi-analytic 
modeling  is  to  adjust  the  model  parameters  to  make  the  computed  electron  density  consistent 
with  the  electric  potential.  The  appropriate  values  of  both  parameters  (ax  and  b)  are  obtained 
by  comparing  the  analytic  solutions  given  in  (19)  with  the  quasi-analytic  solutions  obtained 
by  compressing  the  F-layer  coordinates  using  (42)  for  an  untilted  bubble  or  (67)  for  a  plume 
with  internal  polarization  and  plasma  shear  flow. 

The  vertical  bubble  motion  comes  from  polarization  of  the  horizontal  density  gradients.  The 
horizontal  Pedersen  conductivity  from  (19)  though  the  center  of  the  electric  potential  is  given 
by 


E(‘)(x,0)=  exp- 


a(l  +  m) 


2tan‘ 


B,  -  2  I  x 


S  lb's 


-  sign(a)  7i 


Cp  (l+|x|b)2 
1  -  a  -  B,  |  x  |b  +x2b 


(70) 


where  A,  =^/-a[4b  +  a(b-l)2]  and  B,  =  -a(b  - 1)  -  2  and  |x|  =  Vk^.  Assume  that  the 

electron  density  at  the  equator  is  directly  proportional  to  the  integrated  Pedersen 
conductivity.  In  absolute  coordinates,  the  analytic  model  for  the  horizontal  electron  density 
profile  through  the  center  of  the  bubble  is 


Ne(xp,yc(t))=  exp- 


a(l  +  m) 


2tan‘‘ 


A, 


-  sign(a)  n 


|b\2 


C,  (l+lx/RJ”) 


1-a-B,  |x  /R,|  +|x  /R, 


|2b 


(71) 

where  xp=x-  xc(t)  is  the  horizontal  distance  relative  to  the  center  of  the  potential  function 
and  (xc(t),  yc(t)}  is  the  location  of  this  center.  The  corresponding  electric  potential  is 


„  [x - *.(t)]  ^(t)  +  [y - yc(t)]  VxR (yc(t)) 
®<x'y,t) = B»  _ V  “ 


1+ 


x~xc(t) 

R,(t) 


y-Yc(t) 

R,(t) 


(72) 


where  ax(t),  b(t),  Ri(t),  and  Ci(t)  are  parameters  that  will  be  allowed  to  vary  with  time  as  the 
bubble  evolves.  The  bubble  rise  velocity  V^t)  =  ^Tx  and  the  bubble  retardation 

B0 

velocity  =  ay VxS  have  been  defined  in  the  previous  sections.  The  electric  potential 

follows  a  trajectory  given  by 
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(73) 


-^  =  Vx0(yc)  =  (l-ay)VxS(yc) 

<^_v  ff\  -  ax(t)  Etx 

at  _VyoCt)~  B0 

With  this  trajectory,  the  Lagrangean  coordinate  map 

ax(x0,yo,t)  =  1  a<X>(x,y,t)  d  3yCx0,y0,t)_  1  3Q(x,y,t) 

9t  B0  dy  xS  at  B0  dx 

is  again  used  to  determine  the  distortion  of  the  plasma  layer.  The  normalization  procedure 
fits  the  function  given  by  (71)  to  the  densities  determined  using  the  Lagrangean  transport 
given  by  (74)  as  applied  to  the  stratified  plasma  profile. 

Tests  of  the  normalization  procedure  has  demonstrated  that  the  numerical  value  for  the 
temporal  normalization  constant  ax  ETx/(B0  r0)  is  equal  to  approximately  70%  of  the 

calculated  growth  rate  y  =  (-ETx/B0)/LN  of  the  instability  so  that  the  parameter  ratio 
-r0/ax  »  Ln  /O.7.  The  recommended  procedure  for  providing  reasonable  models  of 

equatorial  bubbles  is  to  first  choose  a  scale  length  ro  that  matches  the  dimension  of  the  “seed” 
needed  to  produce  the  bubble.  Second,  select  the  potential  amplitude  using  the  simple 

expression  ax  0.7  r0/LN  where  LN  =  Ne(y)/ ^Ne(y)  is  the  initial  scale  length  of  the 

dy 

bottomside  of  the  background  ionosphere. 

At  this  point,  all  the  steps  outlined  in  Figure  lb  have  been  completed  and  the  densities  for  the 
ionospheric  bubble  can  be  computed  with  relative  ease.  The  only  numerical  computation  is  a 
solution  of  the  ordinary  differential  equations  (74)  which  are  applied  to  the  model  of  the 
unperturbed  ionosphere.  The  formulas  described  here  provide  the  electron  densities  for 
propagation  of  OTH  radar  signals  through  the  layers  to  determine  the  effects  of  plasma 
bubbles  on  the  ground  clutter  signals.  The  ray  trace  procedure  is  described  in  the  next 
section. 


Vn.  Ray  Trajectories  Through  the  Unstable  Ionosphere. 

Propagation  of  HF  waves  in  the  ionosphere  has  been  studied  for  over  one-half  century.  The 
raytrace  code  developed  to  study  OTH  radar  clutter  is  based  on  the  theory  provided  by 
Hazelgrove  [1954],  Yeh  and  Liu  [1972],  Jones  and  Stephenson  [1975],  and  Budden  [1985]. 
A  ray  trajectory  is  characterized  by  a  position  vector  r  =  (x,  y,z)  and  a  wave  normal  direction 

characterized  by  the  refractive  index  vector n  =  (nx,ny,nz) .  The  wave  normal  vector  k  is 

related  to  the  refractive  index  vector  by  k  =  n  oo/c  where  (0  is  the  wave  frequency.  If  the 
wave  phase  §  along  the  ray  path  starts  with  a  value  of  zero.  At  each  point  on  the  ray,  a 
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dispersion  equation  is  satisfied  with  the  form D(x, y,  z,  nx ,  ny ,  nz ,  co)  =  Constant.  The 
canonical  equations  for  the  ray  path  are  written  as 


dr  VnD  ,  dn  V.D 

—  = - 2 - and  —  = - - - 

dP'  (odD/da>  dP'  (odD/do) 


(75) 


where  Vn  and  Vr  are  the  gradient  operator  in  refractive  index  and  Cartesian  space, 
respectively,  P'  =  ct  is  the  group  path  and  c  is  the  speed  of  light.  The  scalar  multiplier 
dD/dco  is  defined  as 


f®.  =  ^  .  Vn£)  =  —  -  —  •  VnD 

dco  dco  dco  dco  to 


(76) 


The  phase  path  front  P  =  <p —  differs  from  the  surface  of  constant  group  delay  P' .  The  phase 

CD 


front  is  computed  using 


dP  _  dr  =  n  V„D 
dP'~n  ’dP'_  co  dD/dco 


(77) 


Two  forms  of  the  dispersion  equation  are  used  to  trace  the  rays.  The  Appleton-Hartree 
dispersion  formula  is 

D±  =  n2  +  nj  +  n2  - 1  +  X(1  - X)/[l - X - Y2SinV2 ± R(BX , By, Bz,nx, ny  ,n  J]  =  0  (78) 


where  the  radical  R  =  -Jy4  Sin4\|//4  +  Y2  Cos2y  (1  -  X)2  ,  the  plasma  frequency  squared  is 
cd^  =e2Ne/(e0me),  X  =  co^,/co2  the  plasma  frequency  normalized  by  the  wave  frequency, 

Y  =  eB/(meco)  the  gyro  frequency  vector  normalized  by  the  wave  frequency,  and  vp  the  angle 
between  the  wave  normal  and  the  magnetic  field  defined  by  n  •  B  =  n  B  Cos\|/  .  The  +  and  - 
signs  in  (76)  refer  to  the  ordinary  and  extraordinary  modes,  respectively.  The  Booker 
Quartic  dispersion  equation  is 

Dq  =  (1  - X -  Y2) n4  +  X  (n •  Y)2  n2  +[Y2(2-X)-2(1-X)2]n2 -X(n- Y)2  +(1-X)[(1-X)2  -  Y2] 
(79) 

In  the  plasma  where  X  ^  0.1,  dispersion  equation  (79)  is  used  and  if  X  <  0.1,  (78)  is  used  to 
trace  the  rays.  This  selection  preserves  the  identity  of  the  mode  and  allows  the  rays  to 
change  direction  at  reflection  points. 

The  purpose  of  the  ray  tracing  is  to  compute  the  Doppler  shifts  of  the  ground  clutter  signal  in 
a  non-stationary  ionosphere.  Doppler  frequency  shift  along  the  ray  path  is  given  by 
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dAco  _  1  8D/9t  _  1  8D/3Ne  dNe/9t 

dP'  “  c  dD/d©  ~  c  dD/d© 


(80) 


where  it  is  assumed  that  all  temporal  variations  in  the  dispersion  are  due  to  fluctuations  in  the 
electron  density.  For  this  study,  the  electron  density  changes  are  assumed  to  be  the  result  of 
transport  by  electric  fields  as  described  in  the  previous  sections.  From  (1)  and  (22),  the 
electron  density  changes  due  to  incompressible  transport  is  given  by 


=  -v-V,N, 


Vr3>xB 


V  N 

rx  e 


(81) 


Combining  (80)  and  (81)  yields  the  equation  of  the  Doppler  shift  of  the  received  signal  in 
terms  of  the  gradients  in  the  electric  potentials  and  the  electron  densities  in  the  ionosphere. 


dAco  3D/9N  v  _  9D/8Ne  Vrd>xB  _  XT 

- = - — VN  = - - — - - VN 

dP'  dD/d©  c  r  e  dD/d©  cB„ 


(82) 


The  Doppler  shift  for  reflection  from  the  ionosphere  is  referenced  to  the  maximum  Doppler 
shift  for  specular  reflection  from  a  moving  mirror  given  by 


A°>max=2®—  (83) 

c 

where  vy  is  the  velocity  parallel  to  the  line  of  sight.  A  vertical  radio  ray  reflecting  the 
ionosphere  at  normal  incidence  will  have  the  Doppler  shift  of  (83).  For  ground  clutter,  the 
maximum  Doppler  shift  is  2  Aw^  because  for  ground  clutter  the  ray  has  two  reflections 

from  the  ionosphere  and  one  from  the  ground  before  coming  back  to  the  transmission  source. 
A  radio  wave  reflecting  with  oblique  incidence  will  have  a  Doppler  shift  less  than 
2  A©^  with  a  limiting  value  of  zero  shift  for  a  horizontally  propagating  ray  with  no 
ionospheric  interaction. 

The  procedure  for  modeling  the  OTHR  clutter  from  equatorial  bubbles  in  the  ionosphere  is  to 
use  the  quasi-analytic  bubble  derived  in  the  previous  sections  to  describe  the  electron  density 
gradients  and  the  electric  potential  variation.  Rays  are  traced  through  this  ionosphere  using 
(75),  (76),  (77),  and  (82).  The  electron  density  enters  into  the  equations  through  (1)  the 
parameter  X  in  the  dispersion  equations  (78)  and  (79)  and  (2)  the  Doppler  frequency  given 
by  (82).  The  ray  propagates  from  the  ground  radar  source,  is  bent  by  the  bottom  side 
ionosphere,  and  is  terminated  when  the  ray  path  intercepts  the  ground.  The  radar  signal  is 
scattered  by  the  ground  and  the  ray  path  follows  the  identical  trajectory  back  to  the  radar 
source.  The  Doppler  at  the  radar  receiver  is  twice  the  Doppler  at  the  ground  scattering  point. 

As  a  test  of  this  ray  trace  model,  the  Doppler  shifts  in  the  ground  clutter  will  be  computed 
from  the  F  region  vertical  drifts  observed  in  the  quiet  time  ionosphere.  Scherliess  and  Fejer 
[1999]  use  data  from  the  incoherent  scatter  radar  at  Jicamarca,  Peru  and  ion  drift  meter 
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measurements  from  atmospheric  explorer  to  derive  a  global,  empirical  model  of  the  vertical 
drifts  at  the  equator.  These  drifts  are  in  addition  to  the  horizontal  plasma  drifts  illustrated  in 
Figure  14  and  the  internal  plasma  drifts  associated  with  the  equatorial  bubbles.  The 
Scherliess-Fejer  model  provided  the  variation  in  vertical  drift  illustrated  in  Figure  19  for 
calendar  days  100  and  240.  The  largest  variation  in  the  vertical  drift  occurs  between  1800 
and  2000  local  time  near  dusk.  In  the  spring  (Figure  19a)  the  vertical  drift  changes  from  45.7 
m/s  upward  to  37  m/s  downward  over  the  period  of  two  hours.  In  the  summer  (Figure  19b) 
the  variations  in  the  vertical  drift  are  smaller. 
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Figure  19.  Local  time  variation  in  equatorial  vertical  drift  for  quiet  time  periods.  The 
empirical  model  drifts  represent  the  longitude  near  66°  W for  (a)  solar  radio  flux  of 200  for 
a  day  in  April  and  (b)  solar  radio  flux  of  116  at  the  end  of  August. 
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Tracing  rays  though  a  model  ionosphere  Figure  13  with  the  Scherliess-Fejer  drift  model  in 

Q 

(82)  gives  the  predicted  Doppler  shifts.  These  will  be  multiplied  by - ,  giving  the 

'2© 

equivalent  radial  velocity  for  reflection  from  a  horizontal  target.  Figure  20a  shows  the  ray 
paths  which  are  launched  with  zenith  angels  ranging  from  0  to  80  degrees.  The  computed 
Doppler  shifts  in  the  received  echoes  are  illustrated  in  Figure  20b  in  terms  of  equivalent 
radial  velocity.  The  vertical  velocities  for  the  ionosphere  are  taken  from  the  values  in  Figure 
19  at  17:44  local  time.  One  sample  of  measurements  for  the  radial  velocity  at  a  similar  local 
time  is  shown  in  Figure  21.  The  measured  Doppler  shift  of  9.4  m/s  is  shown  by  the  circle  in 
Figure  20b.  The  results  of  this  ray  trace  test  show  that  the  ray  trace  model  yields  reasonable 
results  for  a  uniformly  rising  layer.  In  practice,  the  uniform  lifting  of  the  ionosphere  will 
produce  a  shift  of  all  radar  signals  by  a  constant  amount.  This  shift  will  not  obscure  the 
radial  motion  of  any  targets,  since  they  will  share  the  same  apparent  Doppler  shift. 


One  Way  Group  Path  (km) 


Figure  20.  Model  calculations  of  (a)  ray  trajectories  and  (b)  Doppler  shift  in  ground  clutter 
for  two  vertical  drifts  speeds  at  17:44  local  time.  The  shortest  group  paths  are  for  the 
vertical  rays  with  radial  velocities  twice  the  vertical  speeds.  The  measured  Doppler  shift  is 
illustrated  in  Figure  21. 
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Relative  Power  (dB) 


Figure  21.  Doppler  shift  ground  clutter  on  30  August  2002  at  17:44  Local  Time  in  the 
Puerto  Rican  longitudes.  The  radar  frequency  is  24.4  MHz  and  the  radar  range  bin  is 
around  2500  km.  The  center  of  the  ground  clutter  peak  is  found  at  a  radial  velocity  of  9.4 
m/s. 

Future  research  will  combine  the  equatorial  model  densities  illustrated  by  Figure  18  and  the 
raytrace  code  used  to  generate  Figure  20.  The  equatorial  bubble  densities  will  be  extended 
along  the  magnetic  meridian  using  a  Dipole  field  model  given  by  (51).  With  the  dynamics  of 
the  equatorial  bubbles  used  in  wave  propagation,  complex  Doppler  shifts  will  be  generated 
that  provide  ground  clutter  that  spreads  in  frequency  or  radial  velocity.  The  quasi-analytic 
tools  derived  here  are  uniquely  suitable  for  ray  trace  studies  of  ground  clutter  generation 
because  the  rays  can  be  rapidly  traced  though  the  density  structures  and  because  the  Doppler 
velocities  are  directly  computed  by  projection  of  the  plasma  motion  vector  along  the  ray 
wave  normal. 
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VIII.  Summary 


Using  quasi-analytic  solutions  of  the  plasma  transport  and  radio  propagation  equations,  a 
model  was  derived  to  study  the  effects  of  ionospheric  motion  on  OTH  radar  clutter.  The 
quasi-analytic  formulation  has  been  derived  to  imbed  ionospheric  bubbles  in  a  model 
background  ionosphere.  This  technique  provides  a  three-dimensional  electron  density 
description  using  the  mapping  by  the  electric  potential  functions  of  the  initial  distribution  of 
the  ionosphere.  This  technique  has  the  advantage  over  a  full  numerical  simulation  by 
providing  (1)  factors  of  10  to  100  enhancements  in  computational  speed,  and  (2)  analytic 
formulations  that  are  much  easier  to  us  for  numerical  raytracing.  The  analytic  formulation 
has  been  benchmarked  against  full  three-dimensional  simulations  of  equatorial  bubbles 
[Zalesak,  et  al.,  1982]  currently  in  use  at  the  Naval  Research  Laboratory  for  irregularity 
simulations.  The  effects  on  OTH  radar  ray-path  bending  and  range  Doppler  spread  for 
ground  clutter  have  been  tested  for  vertical  motion  of  the  equatorial  ionosphere.  Future 
studies  will  investigate  for  ray  paths  through  the  equatorial  bubbles  that  are  moving  both 
horizontally  and  vertically.  Finally,  possible  mitigation  strategies  for  reducing  OTH  radar 
clutter  effects  will  be  will  be  tested. 
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Appendices 


Appendix  A.  Numerical  Solution  of  the  Potential  Equation  Using  a  Direct  Solver. 

Numerical  solutions  are  required  when  conditions  of  simplified  geometries  or  uniform  flows 
yield  a  complicated,  nonlinear  partial  differential  equation  for  the  electric  potential.  In 
Cartesian  coordinates,  the  non-separable  elliptic  equation  that  describes  the  electric  potential 
is  given  by  (7)  and  is  repeated  in  equivalent  form  here 


920  92  30  „  920  32  30  92 

Zp  9x2  +_3X"  3X +Zp  39r+_39"  3p"  9x 


(Al) 


The  conductivity  function  2p(x,  y)  is  known  and  the  potential  function  0(x,  y)  is  found  as  a 
numerical  solution  to  the  equation  (Al). 

This  equation  is  converted  into  a  set  of  linear  equations  using  the  usual  finite  difference 
approximations  to  the  derivatives  given  by 


320 

3x2 


Ax2 


320  Q^-2^+0.., 

3y2  Ay2 

30  ^i+1,i-Oi.1,j 
3x  2Ax 

3y  2Ay 


(A2) 


The  Pedersen  conductivity  is  sampled  in  the  uniform  solution  grid  spaced  by  Ax  and  Ay  to 
form  the  array  variables  2pijwith  (i=l,2,  ...,  M)  and  (pi,  2,  ...,  N).  To  complete  the 

solution,  boundary  conditions  of  periodic,  fixed/Dirichlet,  derivative/Neumann  or  mixed 
form  are  provided. 

The  unknowns  Oy  for  (i=l,2,  . . .,  M)  and  (j=l,  2,  . . .,  N)  are  grouped  into  linear  arrays  given 
by 

Xj  =  [Ojj  ...  Omj]  (A3) 

The  resulting  linear  system  can  be  written  as  an  extended  tri-diagonal  matrix 
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(A5) 
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0 

0 

0 
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D2 

D3 

dn_2 

dn-, 

°N. 

XN 

_yn 

where  the  M  x  M  block  matrices  A.,Bj,andCj  are  functions  of  the  finite  difference 

parameters  and  the  Pedersen  conductivity  samples  Zpi  j .  The  matrices  Ai  and  Di  are  needed 

for  periodic  boundaries  in  the  y-direction.  The  string  of  matrices  [Di,  D2,  D3,  . . .,  Dn-i,  Dn] 
allow  inclusion  of  an  additional  condition  on  the  potential  such  as 


J"jd>(x,  y) dx  dy  =  0  . 


(A6) 


This  condition  arises  when  the  addition  of  a  constant  to  a  solution  also  yields  solution.  The 
nonuniqueness  occurs  if  the  boundary  conditions  are  completely  periodic  and/or  specified  by 
constant  derivatives  (i.e.,  Neumann).  With  these  types  of  boundary  conditions,  (A6)  prevents 
the  square  matrix  in  (A5)  from  being  singular  and  a  numerical  solution  can  be  obtained. 
Finally,  the  right  side  of  (Al)  and  boundary  conditions  are  contained  in  the  linear  arrays  Yj . 

The  algorithm  for  solving  (A5)  is  a  generalization  of  the  Thomas  Algorithm  for  scalar 
tridiagonal  systems  [Dahlquist  and  Bjorck,  and  Anderson,  1974].  Initialize  with  new  matrix 
variables 


a,  =  B,1,  S,  =  a,  Y„  a2  =a,  C„  b,  =  -a,  •  Aj  (A7) 

Continue  with  the  equations 

8|  =  Oi_i  •  C|.j  j  ctj  =  [B|  -  Aj  •  a,]  ,  S,=«,  [Y,-A,S,J,  b(  =  —  ctj  •  A(_|  •  b (A8) 

for  the  index  in  the  range  i  =  2, ... ,  N  - 1 .  The  next  operations  define  a  new  set  of  variables 
starting  with  =  0,  b^  =  I ,  where  I  is  the  M  x  M  identity  matrix,  and  continue  with 


^N— I  —  ®N— I  aN-i+l^N-i+l>  ^N— I  ~  ^N-l  aN-i+l^N-l 


1+1 


(A9) 


for  the  index  “i”  in  the  range  i  =  l,...,N-l.  The  solution  for  j  =  N  is  given  by 


XN  = 


-1-1  r 


■>' 


i=l 


•IX  S' 


(A10) 
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The  remaining  solution  vectors  are  found  from 

x,=s;+b;-xN 


(All) 


where  the  index  “i”  is  given  by  the  range  i  =  1, ...  ,N  - 1 .  This  algorithm  was  used  to  provide 
the  numerical  potential  solutions  illustrated  in  Figure  5  and  the  data  given  in  Tables  II  and 

m. 

Appendix  B.  Exact  Solutions  to  the  Elliptic  Potential  Equation  in  Cartesian 
Coordinates 

An  analytic  solution  to  the  non-separable  elliptic  equation  that  relates  the  Pedersen 
conductivity  to  the  electric  potential  in  the  ionosphere  was  given  by  (16).  A  localized  plasma 
structure  may  become  polarized  in  a  gravitational,  neutral  wind  or  electric  field  flow  field. 
The  equation  describing  this  polarization  has  been  given  by  (4)  as 

Vx  •  (EPVX<D)  =  VA  •  J[E0  +  (U+-^-)xB]  apdz  =  ET  (Bl) 


In  Cartesian  coordinates,  (Bl)  was  written  as  (7)  assuming  that  only  a  horizontal  equivalent 
electric  field  Etx  was  present.  Here,  the  total  electric  field  that  drives  the  potential  in  the 
general  case  is  assumed  to  have  the  from 

ET  =  ET(cosa  x  +  sina  y) 


where  x  and  y  are  the  unit  vectors  in  the  x-  and  y-directions  and  a  is  the  angle  between  the 
equivalent  driving  field  and  the  x-axis.  In  Cartesian  coordinates,  this  equation  becomes 


a2o  aLog(Ej  a®  a26  dLog(sp)  ao 


ax2 


ax 


= cosa- 


ax  ay2 

aLog(ip) 

dx 


+  sma 


dy  dy 

aLog(sp) 

dy 


(B2) 


where  the  normalizations  C>  =  <t>/(r0  ET) ,  x  =  x/ r0>  and  y  =  y/r0have  been  used.  A  general 
solution  can  be  found  using  separation  of  variables  with 


0>(X,y)  =  F(X)Gv(y) 


and 


(B3) 


MX,  9)  =  H(y)  Exp[  Lx(X)  ]  (B4) 

With  these  substitutions  into  (B2),  the  solution  for  the  derivative  of  the  Lx  is  found  to  be 
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,  m  _  [  m  Gy(y)  -  sing  ]  H'(y)  +  H(y)  [Gy(y)  F*(x)  +  F(x)  G^(y)] 
xW  H(y)  [cosa  -  Gy(y)  F'(x)] 


(B5) 


The  functions  Gy(y)  and  H(y)are  chosen  so  that  only  a  function  of  x  remains.  Examination 
of  the  denominator  of  (B5)  requires  that 

Gy(y)  =  cosa  (B6) 


to  eliminate  one  of  the  functions  of  y .  With  this  substitution,  (B5)  becomes 


_  tana  H'(y)  -  H(y)  F*(x) 
H(y)  [  F'(x)  - 1] 


The  y  -dependence  in  H(y)  is  eliminated  by  forcing  H'(y)/H(y)  =  m  with  the  substitution 


H(y)  =  Exp(my) 
yielding 


L'x(x)  = 


m  tana  -  F*(x) 
F'(x)-1 


(B8) 


(B9) 


This  result  (B9)  yields  the  second  equation  in  (13)  if  F(x)  =  (a0  +  a,x)  G(|  x  |)  and  a=0.  To 
summarize,  the  one  dimensional  potential  function 

<D(x,  y)  =  F(x)  cos  a  (BIO) 


is  found  from  is  found  from  (B3)  and  (B6)  for  an  external  electric  field  Ex  making  an  angle  a 
with  the  horizontal  axis.  The  corresponding  Pedersen  conductivity  function  is  found  by 
substituting  (B8)  and  the  integral  of  (B9)  into  (B4)  with  the  result 


SP(x,y)  =  erayexp 


rm  tan  a  -  F'(x) 
^  F'(x)  - 1 


dx 


(Bll) 


Rotation  of  the  coordinate  axes  (x,  y)  by  a  so  that  Ex  is  horizontal  will  yield  solutions  for 
obliquely  oriented  potential  functions. 


The  solution  in  (B1 1)  is  only  valid  if  a  *  rc/2  .  For  the  case  that  cos  a  =  0,  (B5)  becomes 


K(V=- 


[  F(x)  Gy(y)  - 1  ]  H'(y)  +  H(y)  [Gy(y)  F^(x)  +  F(x)  G;(y)] 
H(y)  Gy(y)  F'(x) 


(B12) 
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The  only  substitutions  that  eliminate  the  y -variations  in  (B 12)  are 

H(y)  =  1  and  Gy(y)  =  C,  exp(  -  ny)  +  C2  exp(ny) 

where  Ci  and  C2  are  constants.  The  resulting  log-density  function  is 


Lx(x)  =  - 


n2  F(x)  +  F#(x) 
F'(x) 


(B13) 


(B14) 


To  summarize,  for  the  special  case  that  the  Ex  vector  is  aligned  with  the  y-axis  (a  =  tc/2),  a 
potential  function  of  the  form 

<D(x,  y)  =  (C,  e~ny  +  C2  en* )  F(x)  (B 1 5) 

is  produced  by  a  Pedersen  conductivity  variation  with  the  form 


2p(x,y)  =  exp 


rn2  F(x)  +  F*(x)d;, 
J  F'(x) 


(B16) 


Rotation  of  the  solutions  in  (B14)  and  (B15)  by  %I2  yield  solutions  for  a  horizontal  electric 
field  Ej  of  the  form 


<D(x,  y)  =  (C,  e  “  +  C2  e“ )  F(y)  and  IP(x,  y)  =  exp 


rn2  F(y)  +  F*(y) 
J  F'(y) 


(B17) 


These  solutions  co-exist  with  the  horizontal  electric  field  Ej  (a=0)  solutions  of  (BIO) 


0(x,  y)  =  F(x)  and  Zp(x,  y)  =  emy  exp 


f-OiLfc 

JF'(x)-l 


(B18) 


Thus,  with  a  horizontal  external  electric  field  (or  equivalently  a  vertical  gravity  or  neutral 
wind  force),  two  separate  families  exact  solutions  for  potentials  and  densities  are  available 
with  (B17)  and  (B18).  With  rotations  the  equations  (BIO)  and  B(ll),  a  large  number  of 
exact  oblique  solutions  can  be  generated.  The  example  illustrated  in  Figure  2  is  only  a  small 
of  sample  of  the  analytic  results  for  Cartesian  coordinate  geometry. 
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Appendix  C.  Exact  Solutions  to  the  Elliptic  Potential  Equation  in  Cylindrical 
Coordinates 

A  second  general  solution  to  the  potential  equation  can  be  derived  by  using  cylindrical 
coordinates  by  selecting  either  radial  or  axial  symmetric  for  the  function  representations.  As 
defined  previously  in  (4)  and  (Bl),  the  elliptic  equation  describing  the  relation  between  the 
Pedersen  conductivity  (£p)  and  the  electric  potential  (O)  is  given  as 

V±  •  (EpVj®)  =  V±  •  J[E0  +  (U  +  ~)  x  B]  apdz  =  ET  (Cl) 


Assume  that  ambient  forces  yield  motion  in  the  y-direction  by  an  equivalent  electric  field  Etx 
in  the  x-direction  given  by 

Et.  =[Ete  +(U  )B0]  (C2) 

V 

in 

The  cylindrical  (r,  0,  z)  coordinates  are  related  to  Cartesian  coordinates  by 
x  =  r  CosG 

y  =  r  SinO  (C3) 

z  =  z 

and  the  unit  vector  in  the  x-  direction  is 


x  =  cos0  r  -  sin0  0 

In  the  cylindrical  coordinate  system,  the  potential  equation  becomes 


(C4) 


a2o  ao  a2o.  ,  ao  ao 


ar2 


-  +  r— — . 

ar  ao' 


-)+r 


as: 


as„ 


ar  ar  ao  ao 


=  ETx(r2CosO-^-rSin0-f)  (C4) 


3r 


ao 


After  normalization  this  equation  is  written  as 


.2a2o  ^ao  a2o  ^2aiog(i  )36  aiog(E  )ao 

r  — -+  r - 1- — ^  +  r  - - - h- - - - 


a? 


a?  ao2 


a?  a? 


ao  ao 


Aaiog[zp]  _  oaiog[Ep] 

=  r  cos0 - - — rsinG- 


(C5) 


a?  ao 

where  the  previously  defined  normalizations  6  =  d>/ETx  and  r  =  r/r0have  been  used. 


To  obtain  one  class  of  radial-dependent  solutions,  apply  separation  of  variables  with  the 
substitutions 


O(f,0)  =  F(0)Gr(f) 


(C6) 
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and 


Zp(r,9)=H(9)exp[L,(r)]  (C7) 

Then,  solve  for  Lr(?)  in  terms  of  functions  that  only  depend  on  r  .  With  the  substitution  of 
(C6)  and  (C7)  in  to  (C5),  the  solution  for  the  first  derivative  of  Lr(?)  is 

L,  =  [f  sin  e  +  GtfFX  e  )]H/(  e )  +  Gr(f)H(  9  )F'(  8 )  +  f  F(  6  )H(  6 )[  G;(f)  +  f  G;(f)] 

A  ’  ?2H(  0  )[cos  0  -  F(  0 )  G'r(?)] 


(C8) 

The  first  step  in  elimination  of  the  0  dependences  on  the  right  side  of  (C8)  is  to  let 
F(  0 )  =  cos  0  with  the  result 


,  _  Gr(?)  [1  +  tan  0  H'(  0  )/H(  0  )]  -  r  [tan  0  H'(  0  )/H(  0 )  +  G;(f)  +  ?  G;(?)] 

r2[G'r(r)-l] 


(C9) 


Next,  the  0  dependent  relation  tan  0  H'(  0  )/H(  0 )  in  (C9)  is  set  equal  to  a  constant  “m”.  The 
solution  of  tan  0  H'(  0  )/H(  0  )  =  m  is  H(  0  )  =  C,  sinm0  where  Ci  is  a  constant.  With  these 
substitutions,  (C6)  and  (C7)  become 


<!>(?,  0)  =  cos0  Gr (?)  and  Ep (?,  0)  =  sinm0  exp[Lp (?)]  . 
and  (C8)  is  reduced  to 

,  ?  =  r  m  +  (1  +  m)Gr(r)  -  r  G;(r)  -  r2  G;(r) 
r2[l  +  G;(r)] 


(CIO) 


(Cll) 


where  the  right  side  is  only  a  function  of  r .  After  integration  and  substitution  into  (C7),  the 
Pedersen  conductivity  is  found  to  be 


Zp  (r,  0)  =  sinm0  exp  {  f— 


j  m  +  (1  +  m)Gr(?)  -  r  G^)  -  r2  G;(r) 

?2[i+g;(?)] 


(Cl  2) 


Equation  (Cl 2)  identical  to  (16)  if  Gr(?)  =  r  G(r)  and  Lp(r)  =  L' (?)  -  m/? . 

To  obtain  a  complementary  set  of  angular-dependent  solutions,  use  separation  of  variables 
with  the  functions 


O(?,0)  =  Fe(0)G(?) 


(C13) 
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and 


Sp(r,0)  =  H(0)exp[Le(0)] 


(Cl  4) 


and  solve  for  Le(0)  in  terms  of  functions  that  only  depend  on  0  .  With  (C13)  and  (C14), 
(C5)  is  rewritten  as 


L'e(0)  = 


r2  cos  0  H'(r)  +  G(r)  H(r)  F;(00+  r  Fe(00)  G'(?)[H(r)  +  r  H'(r)  ]  +  r  H(r)  G'(7)} 

H(r)[r  sin  9  -G(r)  Fe'(  0  )] 


(Cl  5) 


The  first  step  for  elimination  of  functions  of  r  from  (Cl 5)  is  to  let  G(r)  =  r  so  (Cl 5) 
simplifies  to 


4(0)  = 


r  cos  0 1-T(f)/  H(f)  +  F;(  0 )  +  Fe(  0  )[1  +  r  FT(r)/  H(r)  ]} 
sin  0  -  Fe'(  0 )] 


(Cl  6) 


The  r  dependence  is  removed  from  (Cl 6)  by  letting  H(r)  =  rm  so  r  H'(r)/H(r)  =  m  where 
“m”  is  a  constant.  With  this  substitution,  (Cl 6)  is  further  reduced  to 


4(0)  = 


m  cos  0  +  Fe*(  0  )  +  Fe(  0  )(1  +  m  )} 
sin  0  —  Fj(  0  ) 


In  summary,  if  the  potential  has  the  form 
O(r,0)  =  rFe(0) 

and  the  Pedersen  conductivity  has  the  form 
SP(r,0)=rmexp[Le(0)]. 

This  is  written  out  as  the  integral  equation 


Zp(r,0)  =  rmexp{J- 


m  cos  0  +  (1  +  m)Fe(  0  )  +  Fe'(  0  ) 
sin  0  -  Fe'(  0 ) 


d0} 


(Cl  7) 


(CIS) 


(Cl  9) 


This  solution  is  given  for  completeness.  There  have  been  no  useful  applications  of  (Cl 9)  to 
the  ionosphere  because  the  initial  potential  grows  with  radius  and  the  radial  electric  field 
extends  indefinitely. 


56 


References 


Basu,  B.,  On  the  linear  theory  of  equatorial  plasma  instability:  comparison  of  different 
descriptions,  J.  Geophys.  Res.,  107,  doi:10.1029/2001JA000317, 2002. 

Bernhardt,  P.A.  and  J.U.  Brackbill,  "Solution  of  Elliptic  Equations  Using  Fast  Poisson 
Solvers,"  J.  Comp.  Phys.,  53,  382, 1984. 

Bernhardt,  P.A.,  Cross-B  convection  of  artificial  created,  negative-ion  clouds  and  plasma 
depressions:  low-speed  flows,  J.  Geophys.  Res.,  93,  8696-8704, 1988. 

Bernhardt,  P.A.,  Eye  on  the  Ionosphere:  Ionospheric  Profiling  by  GPS  Receiver  Occultations, 
Taking  Advantage  of  Existing  Earth-Based  Infrastructure,  GPS  Solutions,  9,  174-177, 
2005. 

Budden,  K.G.,  The  Propagation  of  Radio  Waves,  Cambridge  University  Press,  Cambridge, 
1985. 

Chen,  K.  Y.,  H.  C.  Yeh,  S.  Y.  Su,  and  C.  H.  Liu,  Anatomy  of  plasma  structures  in  an 
Equatorial  spread-F  event,  Geophys.  Res.  Lett.,  28,  3107,  2001. 

Huang,  C.  S.,  M.  C.  Kelley,  and  D.  L.  Hysell,  Nonlinear  Rayleigh  Taylor  instabilities, 
atmospheric  gravity  waves,  and  equatorial  spread-F,  J.  Geophys.  Res.,  98,  15,631, 
1993. 

Huang,  C.  Y.,  W.  J.  Burke,  J.  S.  Machuzak,  L.  C.  Gentile,  and  P.  J.  Sultan,  DMSP 
observations  of  equatorial  plasma  bubbles  in  the  topside  ionosphere  near  solar 
maximum,  J.  Geophys.  Res.,  106,  8131, 2001. 

Huba,  J.  D.,  G.  Joyce,  and  J.  A.  Fedder,  Sami2  is  another  model  of  the  ionosphere  (SAMI2): 
A  new  low-latitude  ionosphere  model,  J.  Geophys.  Res.,  105,  23,035,  2000. 

Hysell,  D.  L.,  C.  E.  Seyler,  and  M.  C.  Kelley,  Steepened  structures  in  equatorial  spread-F: 
Theory,  J.  Geophys.  Res.,  99,  8841,  1994. 

Kelley,  M.  C.,  J.  J.  Makela,  B.  Ledvina,  and  P.  M.  Kintner,  Observations  of  equatorial 
spread-F  from  Haleakala,  Hawaii,  J.  Geophys.  Res.,  29,  doi:10.1029/2002GL015509, 
2002. 

Keskinen,  M.  J.,  S.  L.  Ossakow,  S.  Basu,  and  P.  Sultan,  Magnetic  flux  tube  integrated 
evolution  of  equatorial  ionospheric  plasma  bubbles,  J.  Geophys.  Res.,  103,  3957, 
1998. 

Keskinen,  M.  J.,  S.  L.  Ossakow,  and  P.  K.  Chaturvedi,  Preliminary  report  of  numerical 
simulations  of  intermediate  wavelength  collisional  Rayleigh  Taylor  instability  in 
equatorial  spread-F,  J.  Geophys.  Res.,  85,  1775,  1980. 

Kil,  H.,  and  R.  A.  Heelis,  Global  distribution  of  density  irregularities  in  the  equatorial 
ionosphere,  J.  Geophys.  Res.,  103,  407,  1998. 

Otsuka,  Y.,  K.  Shiokawa,  T.  Ogawa,  and  P.  Wilkinson,  Geomagnetic  conjugate  observations 
of  equatorial  airglow  depletions,  Geophys.  Res.  Lett,  29, 
doi:  10. 1029/2002GL01 5347,  2002. 

Ott,  E.,  Theory  of  Rayleigh  -Taylor  bubbles  in  the  equatorial  ionosphere,  J.  Geophys.  Res., 
83,2,066,1978. 

Richtmyer,  R.D.,  and  K.W.  Morton,  Difference  Methods  for  Initial-Value  Problems, 
Interscience,  New  York,  1967. 

Sahai,  Y.,  P.  R.  Fagundes,  and  J.  A.  Bittencourt,  Transequatorial  F-region  ionospheric 
plasma  bubbles:  Solar  cycle  effects,  J.  Atmos.  Terr.  Phys.,  62, 1377,  2000. 


57 


Scannapieco,  A.  J.,  and  S.  L.  Ossakow,  Nonlinear  spread-F,  Geophys.  Res.  Lett.,  3,  451, 
1976. 

Scherliess,  L.,  and  B.  G.  Fejer,  Radar  and  satellite  global  equatorial  F  region  vertical  drift 
model,  J.  Geophys.  Res.,  104,  6829, 1999. 

Sekar,  R.,  R.  Suhasini,  and  R.  Ragvaharo,  Evolution  of  plasma  bubbles  in  the  equatorial  F- 
region  with  different  seeding  conditions,  Geophys.  Res.  Lett.,  22,  885, 1995. 

Stephan,  A.W.,M.  Colerico,  M.  Mendillo,  B.W.  Reinisch,  and  D.  Anderson,  Suppression  of 
equatorial  spread-F  by  sporadic  E,  J.  Geophys.  Res.,  107, 
doi:  1 0.1 029/200 1JA000 162, 2002. 

Sultan,  P.  J.,  Linear  theory  and  modeling  of  the  Rayleigh  -Taylor  instability  leading  to  the 
occurrence  of  equatorial  spread  F,  J.  Geophys.  Res.,  101, 26,875, 1996. 

Weber,  E.  J.,  et  al.,  Equatorial  plasma  depletion  precursor  signatures  and  onset  observed 
south  of  the  magnetic  equator,  J.  Geophys.  Res.,  101,  26,829,  1996. 

Yeh,  H.  C.,  S.  Y.  Su,  and  R.  A.  Heelis,  Storm  time  plasma  irregularities  in  the  pre-dawn 
hours  observed  by  the  low  latitude  ROCS  AT- 1  satellite  at  600  km  altitude,  J. 
Geophys.  Res.,  28,  685,  2001. 

Yeh,  K.C.,  and  C.H.  Liu,  Theory  of  Ionospheric  Waves,  Academic  Press,  New  York,  1972. 

Zalesak,  S.  T.,  Fully  multidimensional  flux  corrected  transport  algorithms  for  fluids,  J. 
Comput.  Phys.,  31, 355,  1979. 

Zalesak,  S.T.,  S.L.  Ossakow,  P.K.  Chaturvedi,  Non-linear  equatorial  spread-F  -  The  effect  of 
neutral  winds  and  background  Pedersen  conductivity,  J.  Geophys.  Res.,  87,  151-166, 
1982. 


58 


